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IMPULSIVE MOTION OF A FLAT PLATE 


By JOHN W. MILES+ 
(Fulbright Lecturer, Auckland University College, New Zealand) 


[Received 8 November 1951] 


SUMMARY 
The force-time history required to produce an abrupt discontinuity in the trans- 
yerse velocity of an infinite flat plate, moving in an inviscid fluid, and the velocity 
and acceleration produced by a suddenly applied force are calculated on the basis 
of the acoustic approximation, using techniques developed for the corresponding 
P diffraction problem. 


1, Introduction 

THE question of transient motion of a rigid body in a compressible fluid has 
been discussed on the basis of the acoustic approximation in a previous 
paper (1) with explicit application to the impulsive motion of a circular 
cylinder. In this paper we shall carry out a similar analysis for a flat 
rectilinear plate. As pointed out in the more extensive discussion of ref. (1), 
the results are directly applicable primarily in the field of acoustics. 

The first problem to be discussed, in which the velocity of the plate is 
prescribed, is formally identical with the determination of the total scatter- 
ing of a transient disturbance by a flat plate, which has been treated by 
Fox (2), or with the determination of the lift distribution on a rectangular 
wing in supersonic flow, which has been treated by Busemann (3), Gunn (4), 
Lagerstrom (5), Stewartson (6), and others, using a wide variety of tech- 
niques. In the following, we shall rely primarily on this antecedent work 
in determining the initial response, but we shall introduce elliptic coordi- 
nates to determine the asymptotic response [for which the analyses of 
tefs. (2)-(5) are inadequate}, following Sieger’s treatment (7) of the plane 
wave diffraction problem. Finally, we shall consider a variational approach, 
which appears to best advantage in attacking the inverse problem of deter- 
mining the velocity response to a prescribed force. 


To render our analysis dimensionless, we choose as a characteristic length 
the half width (/) of the plate specified by y = 0, |x| < 1 and as a charac- 
teristic time (//c), where c is the acoustic velocity. We then specify the 
velocity, q, and the perturbation pressure, p, in terms of a dimensionless 
velocity potential 4, according to 


q(x, ¥,t) = cVd(zx, y, t), 
p(x, y, t) — pc*d,(x, y,t), 
+ On leave from the University of California, Los Angeles. 


{Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 2 (1953)] 
5092,22 K 
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where p is the mass density of the fluid, and a, y, ¢ are all dimensionless: 
suffixes denote partial derivatives. The linearized equation to be satisfied 
by ¢ is then 


rx 1 Py dy- (1.3) 

The boundary conditions to be imposed on ¢ are 
$,(v,0+,t) = v(t) (jz7] <1, t>0), (1.4) 
lim |r'd(r cos 6,7 sin 6,t)| < K. (1.5) 


ae 6 
Equation (1.4) prescribes the velocity of the plate, while (1.5), where K is 
a fixed constant, ensures the appropriate behaviour at infinity. 

In virtue of the symmetry of (1.4) with respect to y = 0, we remark that 
¢ must be antisymmetric thereto, whence, in recognition of the requirement 
of continuity of the solution except across the plate, we may assert that 

f(x, 0, t) Oy (je) > 8). (1.6) 

We next introduce the convention that the Laplace transforms of the 
various time-dependent functions be represented by corresponding capital 
letters; e.g. 


D(x, y, 8) Lid(x, y, t)} [ e-“d(a,y t) dt. (1.7) 
0 


Transforming (1.3), and noting that the initial values of ¢ and 4, vanish 
everywhere except at the plate, where ¢,(7,0+-,0-+-) is finite (vide infra), 
we have the modified Helmholtz equation 
0... | 0, se@D 0 (1.8) 
to which a solution is required subject to the conditions obtained from the 
transforms of (1.4)—(1.6). 
The force that must be applied to the plate to produce the prescribed 
velocity v(t) is given by 


1 
P(t) m(c?/L) f (t) l | [ p(w, 0+, t)—p(a, 0 ,t)| da (1.9a) 
ms 
1 
2pc*l | o(x,0+,t) da, (1.9b) 
nl 


where m is the virtual mass for steady flow, viz. 

m pl. (1.10) 
We seek to exhibit the Laplace transform of the dimensionless force in the 
form 


F(s) = sZ(s)V(s), (1.11) 
where Z(s) is defined by , 
9 ~ 


Z(s)V(s) —— D(a, 0+, 8) dx. (1.12) 
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IMPULSIVE MOTION OF A FLAT PLATE 13] 


Taking the Laplace transform of (1.11) with the aid of the Faltung theorem, 
we obtain t 
f(t) = v(0+)z(t)+ | #(r)z(t—7) dr, (1.13) 
0 
whence z(t) appears as the (dimensionless) force required to produce a unit 


jiscontinuity in the velocity at time zero. 


2. Solutions for small values of the time 

For small values of the time ¢ the solution to the boundary-value problem 
an be described in terms of the disturbances originating at the edges, 
together with their successive reflections at these edges, a technique due 
riginally to Schwarzschild (8), albeit more appropriate to the transient 
problem [see, e.g. Fox (2) and Gunn (4)|. To this end we write 


z(t) > 1(t—2n)z,,(t—2n), (2.1) 
n 0 
where z,(¢) results from the edge waves originating at t = 0, z,(t) from the 


reflection of these waves at the opposite edges, etc. The corresponding 
expression for Z(s) is given by 
Zs) Y e-2"8Z (8). (2.2) 
n=0 
To determine Z,(s) we may consider the reduced problem of a single edge, 
0 e Lo] 
with the modified boundary conditions (after shifting the origin of 2 to the 


“lore 
eage )} 


®,, (2, O+,8) a* (2 > 6), (2.3) 

®, (x, 0, 8) 0 iz <. Wi, (2.4) 
where we have put V’(s) = s~!, following (1.11)—(1.13). The boundary-value 
problem posed by (1.8), subject to (2.3) and (2.4), is essentially the classical 
half-plane problem first solved by Sommerfeld. For our purposes a con- 
venient form of the solution is that given in Lamb’s treatise (9), section 308 
20), after subtracting the incident wave, setting y = 0, and substituting s 


for ik therein; we find 


D,(x, 0+, 8) Fs! | (s€) bef dé, (2.5) 


Z,(8) is now determined by superimposing the two edge effects. Thus, 


setting V(s) s-l in (1.12) and shifting the origin, we have 
1 
Z,(8) wa D(x, 0-+,s8) dz. (2.6) 
9 


Sube**Suting (2.5), taking the inverse transform, and integrating with 
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respect to 2, subject to the temporary restriction t < 
effectively with only half the plate), we obtain 


1 (since we are dealing 


Z(t) —(1—}#). (2.7) 


The determination of z,(f—2) follows by analogy from the work of Gunn 
(4) or Lagerstrom (5), whence we obtain 








ba, | (#2 + 
2,(t—2) = = (Ca -1)cos-1(2/t) + 4(2—4)! nf 2", (2.8) 
the usual principal values being understood. 

FZ (t) 
1:0 
0:8 
0-6 
0-4 
0-2 

—iuaiiiialiitas: t 


Fic. 1. (47)z(t) (2.1) and (4.9)], the (dimensionless) force 
required to produce a unit step in the velocity. 


[as determined by 


fa bl . . ° oe, fe ° . 
Che subsequent expressions for z, are unwieldy, and it is simpler to seek 


n 


an asymptotic solution for ¢ > 4. 
7) and (2 


The results so era [ (4.9) infra} 


together with those of (2 8), are presented in Fig r 


3. Solution in elliptic coordinates 


An alternative approach to the boundary-value problem is afforded by 


the introduction of the elliptic coordinates | 

x = cosh é cos », y = sinh ésin yn, (3.1) 
the strip |x| < 1 being specified by € = 0. In these coordifates ( .8) trans- 
a P ‘ { 
forms to O.-+ ®,, —s*(cosh*é—cos*n)® = | (3.2) 
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IMPULSIVE MOTION OF A FLAT PLATE 


while the boundary conditions (1.4) and (1.6) transform to 
®; V(s)sinyn (€ 0), (3.3) 
® 0 (y 0 and z=). (3.4) 


Asolution of (3.2) satisfying (3.3) and (3.4), together with the appropriate 


onditions at infinity, is given by 


= V(s) ¥ (—)"AP"M(g)[Gek,,, .4(€, —g)/Geky, 410, —9)}8€ on 41(, —9): 
n 0 

(3.5) 

q 8/4, (3.6) 

where se, and Gek,, are those solutions of Mathieu’s equation analogous to 

inn@ and K,,(sr)in the polar coordinates r, 0, and the A‘)(q) are the Fourier 

efficients in the expansion of se,,(7, —q) in sinmy. This type of solution 

sdue essentially to Sieger (7), although our notation is that of McLachlan’s 

treatise (10) (hereafter denoted by ) followed by the appropriate section 
number. 

Setting € — 0 in (3.5), substituting in (1.12) and integrating, we obtain 


Z(s) = > Fau(Q[APr@s, (3.7) 
F. (q) Gek,, (0, —q)/Gek’,(0, —q), (3.8) 


where accents denote derivatives. 

In principle, this result is valid for all s, but the convergence of the series 
isslow for large values of s. However, in the solution of the original problem 
for large values of ¢ it suffices to determine the Laplace transform for small s, 
while for small ¢ the solution given in the preceding section is available. 

For the subsequent development of F,, in powers of q it is convenient to 
ntroduce the representation | M 13.30 (7)| 


(se) -+J,,,(38e-£).K,(38e8) 
(3.9) 


Gek . 


i1(§, —) C, > APA @LLGse ‘)K 


r=U 


r 


where J, and K, are the modified solutions to Bessel’s equation, and C,, is 
i constant, whose value is not required. Differentiating (3.9), making use 
if the Wronskian for (J., K.) and substituting in (3.8), we obtain 


Pana) = ( ¥ AMP (@){(2r+-1)4+ 2¢[T, (48) K-48) — J, 4148) K7 (38) ]} ) 


2 ARMs (q). (3.10) 
4. Asymptotic solution — 

The coefficients A{?”*)(q) can be developed in power series with a radius 
ol convergence that is probably at least 4 (cf. M 2.17) but appears 
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empirically to be closer to 4. a required expansions may be determined 





with the aid of M 2.151, 2.21, 3.10, and 3.33, and we abtain the following 
results: . ia vor i 
[A{?(q)}? = 1—(¢/8)?—2(¢/8)3+ Of (q/8)*]. (4.1 
[A{(q)]? = (¢/8)?+2(q/8)3+ Of (q/8)4]. (4.2) 
: l 8)—8(q/8)?+ 4 O} (¢/8): 
F(q) = | —(q/ amt r+) ger. [(q/syA] as 
(1+ (5—8p Yq /8)-4 (3—S8r)(q 8)? (73 ig — 3 ¥)( 8)? + O| (q¢/8)*v]} 
: 1+-(n?+-n)-1(q/8)- a (q/8 21 | 
F, ) 2n+-1 1 HE aria SPR) es Bits ( my 
anil : 1-4 3(n2-9 n)(q/8 ma = a ), (44 
Vv y+In(s/4), (4.5) 
where y is Euler's constant (0-577...). 
It is to be assumed that each of the F,,,,, has at leas it one pole, so that 


Z(s) has an infinite number of them, but the expansions of (4.3) and (4.4) 
serve only to determine the smallest of these with acce Pri able accuracy, viz. 
[cf. Stewartson (6) | 


Gek;(0, —q)|,-s.5 = 9; ™ 0-682+-7 1-114 + 1-3le%?129, (4.6) 


ar 


where the bar denotes the complex conjugate. On separating from Z(s) the 


singularities associated with these poles and neglectirig those of higher 
order, we may expand in ascending powers of s, with the final result 


4 i 4 
Z(s) a_,|(8—8p) 1 es, * z (s 9)” | Lay i] (sa) 14 So 1 > (s 8) 

0 m 0 
1+ (v—$)(s/2)? + (v?—{v-+ y6)(8/2)*4+- O(s8 In’s), (4.7) 
a_, = 0-760exp(—i0-723), | (4.8) 
(The terms of order s* were included in the calculation df the residue, @_,, 


but not otherwise retained.) 
The inversion of (4.7) furnishes the asymptotic expansion 
2(t) 2 re(a_, e%’)—ht-84 3{In(4f) 33 \t 54 Olt 7 |n22) (4.9) 
to which terms like s” make no contribution. [The inve rsion of terms like 
m(y+-Ins)” is discussed in ref. (1).| 
While the positive powers of s in Z(s) make no a a to 2(t) for 
t > 0, we remark that they will contribute to the inverse ‘transform of sZV 
if V contains terms like s-”, m > 1. Thus, if we seek the force required to 
produce unit acceleration, so that V = s~?, the leading term in the expan- 
sion of sZV is s~!, corresponding to a leading term of unity in the inverse 
transform of sZV for large t. We also remark that, while in principle we 
are not justified in retaining the exponential terms and at the same time 
neglecting terms like ¢~’ In*t, numerically the exponential, term is found to 
be dominant throughout much of the useful range of the approximation. 
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IMPULSIVE MOTION OF A FLAT PLATE 135 
Within the (slide rule) accuracy of the numerical calculations, (4.9) agrees 
with the result given by (2.1) over the entire interval 3 < t < 4 and there- 


fore furnishes the desired continuation for t > 4. 


5. Inverse problem 
We consider next the inverse (and physically more interesting) problem 
of determining the motion produced by a prescribed force. If y(t) represents 


the velocity response to a unit discontinuity in the force at time zero, then 


t 
v(t) f(O0+-)y(t)4 | f(r)y(t rt) dr, (5.1) 
0 


ind comparing with (1.11) we have 


Y (s) s—?Z-1(s). (5.2) 

Considering first the case of small values of ¢, we write, as in (2.1) and 
y(t) > 1(t—2n)y,,(t—2n), (5.3) 

n 0 
) S) > e ney’ (a). (5.4) 
n—0 
Substituting (2.2) in (5.2) and expanding in, and equating powers of, e-*°, we 
obtain . On — 
ul Y,(s) s-2Z 1g), (5.5) 
Y,(s) s~*Z5 *(s)Z, (8). (5.6) 
We remark that Z,(s) must be obtained from (2.7) [rather than (2.6), since 
latter represents a problem that is truly equivalent only for ¢ < 1], 
whence 

4 io 
Z (8) s—1(1—}s-1). (5.7) 


Substituting in (5.5) and inverting, we obtain 
Yo(t) Le", (5.8) 


Turning to y,, the substitution of (5.7) in (5.6) yields 


Yi(s (47)?(s— })-*s?Z,(s), (5.9) 
the inverse transform of which is [since z,(0) = 2,(0) = 0] 
t 
y,(t—2) (far)? | (t—r)e#-,(7) dr (5.10 a) 
t 
‘ 7+2\} 
} T + MC r)e'!-) dr, (5.10 b) 


where Z, has been evaluated from (2.8). 
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The integral in (5.10) is intractable, but it may be expanded in ascending 
powers of (¢—2)!, the first term being given by 


Y(t 2) —i(t- 2)3 + O[ (t—2)], (5.11) 


whence it appears that the curvature of y(t) is discontinuous at t — 2 


To determine y(t) for large ¢, we substitute (4.7) in (5.2) to obtain the 


expansion Y(s) s2—1(v—3)—A(v })s?-4+ O(sty?). (5.19) 






1:0 





i i... 








0 





+ ¢ 68 5 6, 7 8 


Fic. 2. y(t)—t and y(t) [as determined by (5.8) and (5.13)], the difference between the 
(dimensionless) velocity and its asymptotic value and the (dimensionless) accelera- 
tion respectively, due to a unit step in the applied force. 


| While we have neglected possible poles in Y(s), the numerical results imply 
that, if present, they are unimportant to our discussion, in contrast to those 
in Z(s).] Inverting (5.12) yields the asymptotic result 

y(t) t+}t-14 Jt-34 O(t- Int). (5.13) 

Att = 1 and 2, (5.13) yields 1-289 and 2-129, respectively, for the values 
of y(t) compared with the exact values 1-294 and 2-135 given by (5.7). Hence, 
in so far as y(t) [but not ¥(t); vide infra] is concerned, it appears that the 
asymptotic result is satisfactory for ¢ > 2, albeit failing to predict the 
singularity at ¢ = 2. 

The numerical results, based on (5.8) and (5.13), are plotted in Fig. 2. 
The curve y(t)—?t represents the difference between the (dimensionless) 
velocity and its asymptotic value, while #(t) represents the acceleration 
(vide infra, section 7). 
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6, Variational formulation 

A rather different approach to the boundary-value problem formulated 
in section 1 is afforded by Schwinger’s variational method (11, 12). 

We proceed by writing the solution to (1.8) in the form 


1 
O(x, y, 8) - | dr exp| iA(~—£&) — (A? +8”)! |y | |@(E, 0+, 8) d&, 
- “1 (6.1) 


where the phase of the radical is chosen to ensure a positive real part, 
thereby satisfying the boundary conditions of (1.5). Equation (6.1) reduces 
tothe required identity at y = 0+ (the negative sign being taken for y < 0) 
in virtue of Fourier’s (integral) theorem. To determine ®(x,0+-,8) we 
invoke the boundary condition (1.4), which yields the integral equation 
1 
V"(s) | G(a—€,s)M(E€,0+-,s8) dé (\x| < 1), (6.2) 


“4 


G(x, 8) lim ; (A? +-s®)! exp|iAw— (A?+-8?)!y| dA. (6.3) 


y—>0 asi 


To cast Z(s) in variational form, we multiply both sides of (6.2) by ® dz, 


integrate over |a 1, and divide by (1.12), whence 
9) 1 2; 1 1 - ' 
Z(s) ca ® dx / | | OG® dirdé. (6.4) 
7 | 4 


It may be shown that (6.4) is stationary with respect to first-order variations 
of ® about the true solution to (6.2). Moreover, if we take the reciprocal 
if (6.4) and substitute (6.3), we may write 
- : \2 
Z-\(s) + lim | (A?+-s?)!exp[—(A?+-s?)ty]| | DeiAx de/ | Mdx\ dy, 
0+ - d ‘ 
1 “1 
(6.5) 
which establishes the positive definite character of Z(s) for real s. (For 
complex s, the real part of Z~! is positive definite, while the imaginary 
part has the same sign as the imaginary part of s.) Finally we remark that 
6.4) is invariant with respect to a scale transformation of ®. In view of 
these properties, (6.4) furnishes an ideal vehicle for an approximate solution 
ot our problem. 
Returning to the kernel, G, carrying out the integration in (6.3) yields 


Res - ~ . . 
G(x, 8) lim h. ; K,[ s(x?+-y?)!], (6.6) 
y—>o+ Oy? 


>0 


7 
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where the positive phase of the radical is chosen. Separating out the 
singularity at the origin, we write 


G(x) G(a, 0) lim a-!(a?+-y?)-?(y?—2?), (6.7 
y—0 
: . ] " " 
G(ax.s) G(x) +—|x-*—s|x|-1K,(s\x])] (6.84 
; to 
Ga) + | K,(o'x)|)o do, (6.8 b 
0 


where the subscript a implies ‘asymptotic’, in recognition of the fact that 
G, determines the behaviour of the result for large ¢. The corresponding 
form for Y, as given by (5.1) and (6.4), is 


8 1 1 1 ' 
Y(s) = ¥,(s)- 2. | odo | | O(s, x) Ky(o|a—E|)@(s, €) dade | | ® dx'\ . 
s. J Ud 
(6.9 
-_ . 2 
Y,(8) = (7/28) | | O(s,x)G,(x—€)O(s, €) das | DM dx|. (6.10) 
1-1 ‘—1 ' 


In order to illustrate the utility of the variational formulation, we 
approximate ® as constant across the plate, which it is exactly only att = 0. 
Actually this approximation is unsatisfactory for Y, (which, however, makes 
no contribution at ¢ = 0) due to the order of the singularity in @,, but it 
may be established directly (from the concept of virtual mass) that y, = t 
| ef. (5.16)], in so far as the dependence of ® on s and z is separable. Then, 


using the approximation only in the second term of (6.9) and introducing 


the change of variable » = o ~—&\, we find 
Y(s) Y_(s) | (252)-1 ( do ( (1 lo In) Ko(n) dy (6.lla 


0 0 
Y(s)+8-! | Ko(20) do+ (2s)-*{ 28K,(2s)—1—K,(2s)—y—Ins]. 
0 (6.11 b) 
Inverting (6.11 b) and taking y, = ¢, we find 
y(t) to +4i(2+Int) (t < 2) (6.12 a) 
4 sin-1(2/t)+-4[t+ (#2—4)!]— 1¢tIn[$4+-34-(2—4)!]  (t > 2). 
(6.12 b) 
In the limit when ¢ (0, (6.12 a) is in agreement with (5.7), as it must be 
in consequence of the approximation used; however, the terms of order f 
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are not in agreement because the value of y,(t) has been chosen to agree 
with the known result for large ¢ rather than small /. (In general, a con- 
sistent variational approximation, i.e. one that poses the same ® in both 
terms of (6.9), would give the exact coefficient of s*” if the value of ® chosen 
were correct to order s?”.) For large ¢, expanding (6.12 b) yields 


y(t) = t+ }t7+3t7+O0(t-), (6.13) 


which is to be compared with the result (5.13). 

Fort 2, (6.12) gives y 2-132 compared with the exact result of 2-135. 
For larger values of ¢ we may expect the agreement to be even better, since 
the first two terms in (6.13) are exact, while the term in ¢-* is not too wide 


of the mark 


7, Acceleration under constant force 
Ifa constant force P, is applied to the plate at ¢ = 0 the resulting (dimen- 


sional) acceleration is given by, ef. (1.9) and (5.1), 


a(t) = (P,/m)y(t). (7.1) 

f we express the apparent mass as a function of time, then 
a(t) P,/m(t), (7.2) 
m(t) m/y(t), (7.3) 


where m(¢) appears as the ‘dynamic virtual mass’. 


The acceleration function, as determined by (5.8) and (5.13), is found to be 


y(t) lre!’ (0 <t < 2), (7.4) 

y(t) 1—}¢-*— St-44 O(t-* Int). (7.5) 

Equations (7.4) and (7.5) are sufficient to determine 7(¢) except for small 

positive values of ((—2). In the latter neighbourhood we may expand #,(¢) 
ibout # 2 and utilize the approximation (5.11) to obtain the result 

y(t) we/16)t—4(t—2)!+ Of (t—2)?]. (7.6) 

It appears that the interaction of the edge disturbances for t = 2 produces 


in asymmetric cusp in the acceleration, a phenomenon that has no parallel 
in the case of the circular cylinder (1). Further, the dynamic virtual 
mass m(t), while exceeding the static value (m) over most of the range of ¢ 
by a factor of 8/7 att — 0), actually falls below this value in the neighbour- 
hood of ¢ = 


REFERENCES 
1. J. W. Mites, Quart. J. Mech. and Applied Math. 4 (1951), 388. 
2. E. N. Fox, Phil. Trans. Roy. Soc. A, 241 (1948), 7. 
3. A. BuseMANN, Jahrb. de? Luftfahrtforschung, 7B (1943), 105. 



























IMPULSIVE MOTION OF A FLAT PLATE 


4. J. C. Gunn, Phil. Trans. Roy. Soc. A, 240 (1947), 327. 


5. P. A. LAGerstTroM, Douglas Aircraft Co. Rep. SM—13110 (Santa Monica, Calif 
1947). "| @O 
6. K. Stewartson, Proc. Camb. Phil. Soc. 46 (1950), 307. 
7. B. Srecer, Ann. d. Physik, 27 (1908), 626. 
8. K. ScHwarzscHILD, Math. Ann. 55 (1902), 177. 
9, H. Lams, Hydrodynamics (Dover, New York, 1945). By 





10. N. W. McLacuian, Theory and Application of Mathieu Functions (Oxford, 1947), 
11. H. Levine and J. Scuwincer, Phys. Rev. 74 (1948), 958. 
12. J. W. Mites, J. Appl. Phys. 20 (1949), 760. 


cons 
pat 
not 
Ex 
circ 
to 1 


the 











ON THE SLOW MOTION OF AN ELLIPSOID IN A 
ROTATING FLUID 
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First received 20 December 1951; resubmitted 21 February 1952] 
Part I. A ForMAL SOLUTION 


SUMMARY 
In the first part of this paper we consider the perturbation in the velocity of an 
iscid 
w motion, relative to the fluid, of an ellipsoid having one of its axes parallel to 


ljuid, rotating about an axis with angular velocity 4, which arises from the 


ixis of rotation Oz of the fluid. Formulae in the form of Bromwich integrals 
found, giving the perturbation at any time, and in the second part we use them 
find the ultimate velocity distribution when the ellipsoid is moving with uniform 
ity. This velocity is steady in general and markedly different according as we 
nsider points inside or outside the cylinder C, which circumscribes the ellipsoid 
nd has its generators parallel to Oz. If the ellipsoid moves parallel to Oz, then 
nside C the fluid is pushed along in front of the sphere as though solid, in agreement 
experiment, while outside C there is a shearing motion parallel to Oz. If the 
psoid moves in a plane at right angles to Oz, then inside C the stream-lines lie on 
ngruent ellipsoids, while outside C they lie on planes perpendicular to Oz but their 
ttern there is rather complicated (see Fig. 2). The agreement with experiment is 
t so good as before, but an explanation of the discrepancy is easily formulated. 
Exceptions to this general picture occur on C, on the ellipsoid, and, in certain 
reumstances, on the axis of C. Finally on C the velocity becomes infinite according 
the linearized theory, and in Appendix C we indicate how viscosity would modify 


flow there 


1. Introduction 
We are concerned here with the slow motion after an impulsive start of an 
ellipsoid in a rotating fluid when one of its principal axes is parallel to the 
axis of rotation of the fluid. It is supposed that the ellipsoid is initially at 
rest relative to axes rotating with the fluid, and that at time ¢ = 0 it is 
given, impulsively, a velocity of translation and that thereafter it moves 
with prescribed velocity and without any relative rotation. The general 
solution is found in this part, and in the second part we use it to deduce 
the ultimate velocity and pressure distribution when the ellipsoid moves 
with uniform velocity both along and at right angles to the axis of rotation. 
The characteristics of the velocity distribution caused by the motion of 
a finite body through a frictionless fluid are greatly different according as 





the fluid is at rest at infinity or is rotating about an axis there. Thus if the 
body is a sphere and the fluid is at rest at infinity, the disturbance caused 


[Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 2 (1953)] 
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by the sphere dies away smoothly as we move away from it. Further, the 
flow is everywhere irrotational and dependent only on the instantaneous 
velocity of the sphere, being steady if the sphere is moving with uniform 
velocity. Ifthe fluid is rotating rapidly about an axis, however, the situation 
is entirely different. The perturbation in the fluid velocity depends not only 
on the instantaneous velocity of the sphere but also on its past history, and 
in general is neither steady nor irrotational. The ultimate motion of the fluid 
when the sphere is moving with uniform velocity is particularly striking. 

If the sphere is moving parallel to the axis of rotation, the cylinder (, 
circumscribing the sphere and having its axis parallel to the axis of rotation, 
separates out regions with markedly different velocity distributions. Pro- 
fessor Taylor (1) found, experimentally, that the fluid inside C appeared to 
be pushed along in front of the sphere as though it were solid. The author (2) 
has confirmed this observation theoretically and has found that in addition 
there is a swirling motion, singular on C’, about the axis of C. Further he 
has found that outside C there is a shearing motion parallel to the axis of 
rotation, singular on C and tending to zero as the distance from C tends 
to infinity. The motion is two-dimensional, i.e. the same on all planes 
perpendicular to the axis, and the sphere experiences a resistance to its 
motion. Exceptions to this theoretical picture occur on the axis of rotation 
and on the sphere, where the fluid oscillates, and on C where it increases 
indefinitely. 

Taylor (3) has also observed the effect of moving the sphere in a direction 
at right angles to the axis of rotation. He found that the fluid inside C 
appeared to move with the sphere while C behaved like a solid obstacle to 
the fluid outside it. Theoretical work on this motion has been carried out 
by Grace (4) using expansions in powers of ¢t. His calculations suggested 
that the ultimate motion was neither steady nor small in the neighbourhood 
of the sphere. From the calculations of this paper it is possible to deduce 
all Grace’s results, and we conclude that his computations were correct but 
misleading. The ultimate flow as predicted by the assumption of small 
motion will be discussed below. 

Taylor (5) has also discussed the motion of a cylinder through a rotating 
fluid, and has found that the stream function of the relative motion may 
be deduced from the appropriate stream function when the fluid is not 
rotating and that the fluid behaves in a similar way. Further, the motion 
only depends on the instantaneous position and velocity of the fluid. 

The purpose of this paper is threefold. Firstly to obtain the ultimate 
fiow pattern as predicted by the assumption of small relative motion when 
a sphere moves perpendicular to the axis of rotation, and to compare this 
with Taylor’s work. The main features of this theoretical solution are that 
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inside C the relative stream-lines lie on spheres with the same radius as the 
moving sphere and with centres on the axis of C, while outside C the stream- 
ines lie in planes perpendicular to the axis of rotation, and their pattern 
there is rather complicated. The agreement with Taylor’s observations is 
not so good as before, but at least, if due allowance is made for the difficulty 
xperienced by real fluids in turning sharp corners, it is possible to offer an 
xplanation of his results on the basis of this solution. Exceptions to the 
theoretical solution sketched out above occur in the same way as when the 
sphere moved parallel to the axis of rotation. There is also a force on the 
sphere inclined to its direction of motion by about 142°, and the motion is 
the same on all planes perpendicular to the axis of rotation. 

Secondly we wish to determine whether the characteristics of the motion 
jue to the sphere are reproduced with more general bodies, and finally we 
wish to bridge the gap between the complicated motion caused by the sphere 
ind the simple motion caused by a cylinder. For these purposes we consider 
nellipsoid with axes of length 2a, 2b, and 2c, of which one axis, of length 2c, 
is parallel to the axis of rotation of the fluid. It is possible to consider a 
arge variety of ellipsoids by varying 6 and ¢ relative to a, from the elliptic 
ylinder to the disk and from the sphere to the thin rod of finite length. We 


| shall find that if the ultimate motion of the ellipsoid is steady the ultimate 


field of flow is also steady, in general, and very like that already described 
for the sphere. The cylinder circumscribing the ellipsoid and having its 
generators parallel to the axis of rotation, which we shall also call C, plays 
the same important part in the ultimate motion. The main difference, in 
fact, occurs on the axis of C’, where the ultimate motion is only oscillatory 
fa =b. The way in which the ultimate motion depends on c is also of 
nterest. If the ellipsoid is moving parallel to the axis of rotation then the 
iltimate flow is independent of c, so that even if the ellipsoid is very 
elongated the ultimate flow is the same as that for a disk; clearly, however, 
fa cylinder moves parallel to itself there should be no effect on the fluid 
lue to its motion. If the ellipsoid is moving at right angles to the axis 
{ rotation, the ultimate flow does depend on c but only in a relatively 
simple way. As coo the pressure on the ellipsoid tends to the two- 
limensional value, but the ultimate flow patterns are not identical, and in 
iact the limit of the three-dimensional velocity distribution has a singularity 
n C and is rotational. We can show, however, that as c increases, the 
insteady terms in the solution oscillate with frequencies which cluster more 
nd more closely round zero, and this implies that it takes longer to reach 
sensibly steady ultimate flow. 

In Appendix C a erude estimation of the modification caused by viscosity 
in the neighbourhood of C is made. According to this estimate the velocities 
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on C are not singular but have maxima of the order of V R!, where V jg q 
representative velocity and R a Reynolds number. 


2. The equations of motion 

We shall take the angular velocity of the fluid to be $, and a rectangular 
cartesian system of coordinates in which the origin O is the centre of the 
ellipsoid, Oz is parallel to the axis of rotation and along one of the principal 
axes of the ellipsoid, and Ox, Oy are along the other two principal axes and 
rotate about Oz with the same angular velocity as the fluid. 

Relative to a parallel set of axes O'x’y’z’ having its origin fixed on the 
axis of rotation, let the coordinates of the centre of the ellipsoid be (z’, y’, z’), 
its velocity components be (U, V, W), and the components of fluid velocity 
be (u,v,w). Further, let the equation of the ellipsoid be 


A x a2 a y? 2 + 22 c2 = ] (2.1) 


relative to Oxyz, where, for definiteness, we shall take a > b. Provided 
that the velocity of the ellipsoid is sufficiently small for squares and products 
to be neglected in the differential equations of motion, these may be 
linearized and we get 


aD 7 “me 
ou, oF 0. a _9 
ot Ox ot | cy 
: (2.2) 
ew . oP Cu. ov. ew 
—— of . —O mé = ies cone ena 
7 2 Ox Oy @ 
where, if p is the pressure and p the density at any point of the fluid, 
p = pP+4pf{(e+a’)?+(y+y’)t. (2.3) 


We cannot justify this linearization of the equations of motion in- 
definitely, because they predict that the vorticity ultimately oscillates 
infinitely, but it is believed that the quadratic terms do not have a significant 
effect on the ultimate flow except possibly near the singularities. 

In the problem under consideration the ellipsoid is given, at time t = 0, 
an impulsive velocity and thereafter moves with a velocity whose compo- 
nents are (U,V, W) relative to O’x'y’z’ so that its principal axes always 
remain parallel to O’x'y’z'. The boundary conditions are that 


zu yo zw x2U yy 2W 
— % +— =-—,+ Y 4° when A =1,t>0; 
az bb? @ a” b2 c2 
u »y=w=0 whnA>l.i 0; 
and u,v,w>0as A>ow,t>0, (2.4) 


and express the condition that on the ellipsoid the component of the fluid 
velocity normal to the ellipsoid should be equal to the component of the 
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velocity of the ellipsoid in the same direction; that before the ellipsoid starts 
tomove the fluid is at relative rest, and that at any time the fluid infinitely 
far from the ellipsoid must be at rest relative to the coordinate axes. 


Let P(x, ¥, 2,8) ( e-“ P(x, y,z,t) dt, (2.5) 
0 
wy ae 
so that P(x, y, 2, t) —— e“P(x,y,2z,8)ds (y> 0). (2.6) 
aT 


We may define a, ¢, @, U, V, and W in a similar way, and may then 
restate the problem as follows. We have to find a function P satisfying 
er or e?P 1-4 s? 


ox* C y* Cz* s* 


0 (: 


bo 
~I 


and such that 


S| : = 
\a* C2 b- cy 


[x oP Ape OE kelp ¥ \+2 a cP 


a= ¢ y b? Cx 
of > Wy 
(924 1)(" — y! 7) when A = 1 


cP oP cP 
, and >O as A> oO. (2.8) 
Cx cy C2 


Once P has been determined we may deduce @, @, and #, for 


s @¢éP l oP 


i 2.9 
° etléx s*+1 oy \ 
D . “ap 
‘ Ll. @F __ 8 ol , (2.10) 
s-t+] or s?4 l cy 
in it pall (2.11) 
S C2 


3. The formal solution 

It is clear that if we replace z by 82z(s? +-]) ‘and c¢ by se(s?+ 1) * the 
problem is then of a type similar to that of the motion of an ellipsoid in an 
necompressible fluid at rest at infinity and it may be formally solved in the 


same way (6). Consider 


—-4 onan ~ 3.1 
6+a? 6+ h2 we J s?) s?] | 2 ( ) 


which defines a pseudo-confocal system of quadrics. Through any point 
wa.28 L 
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(x,y,z) there pass three quadrics characterized by 6 = A, @ = p, and @ 


Vv. 


Hence we may express 2, y, z in terms of A, », v obtaining 
(1+-s?)(A+a?)(u+a?)(v+a?) 

(a? b*)| (1-4 s?)a? s*c?] 
(1 +-s?)(A+-b?)(u+5?)(v+-6?) 

(b2 a*)| (1- 8*)b? s?c?| 


and Z 


| 
9 [(] } s?)X { s’c?][ (1 { s*)u + se? + 8*)p | 8*c?] 
a>( 1-4 8?) s*c? || b> ]+ s?)— s*c?| 


9 
8° 





Let us take A = 0 on the ellipsoid A 1 and suppose that if s is real A + & 
with x, y, and z. Then A, which satisfies 


F(x, y,2z,A, 8) 
[ (1+-s?)A+ s*e?](A+a?)(A+ 6?) — 2°[ (1 +-s?)A+ °c? |(A-+- 6?) 
y?| (1+-s?)A+ s?e?|(A+- a?) — s222(A+a?)(A+6?) = 0, (3.3) 


is defined as a single-valued function of 2, y, z, and s if we insert cuts in 
the s-plane at the branch points of A qua function of s, since and v are both 
negative when s is real. The branch points are shown in Appendix A at the 
end of Part I of this paper to lie on the imaginary axis of s between the 
points s 7. In Part II we shall need the value of A when s is small, and 
from equation (3.3) we see that when s = 0 


(A+ a?)(A+b?)—2?A\(A+ b?)— y?2A(A+a?) = 0, (3.4) 
so that when s 0 either A 0 or 


2A = x?+ y?—a*—b?+ | (x?+ y?—a?—b?)?+ 4(22b?+ y2a?—a*b?)|! = 2 


since A is positive when s is real. 
Now A, vanishes only when 2x*b?+-y?a?—ab? = 0, i.e. at points of the 
circumscribing cylinder C, and 


ne AR (3.6) 
when x = y = 0; further 
(ett y2M tat) tater 


Am Ss (3.7) 
]-+-s? 

when x, y, and z are large. Hence, since A is a continuous function of s 

when s = 0,A =A, if 2b?+-y2a? > a®b? and A = 0 if xb? +y?a? < ab’. 

Inside C, therefore, when s is small 


9.9 
s<2e 


A — . — se2, (3.8) 
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In terms of A, uw, and v the differential equation satisfied by P is 


(pA v)) (a" 1 A)*(b? \)4] se? +- (1-4 s*)A|! . \"P sec oe. = VY, (3.9) 


eA} 
‘he other two terms being obtained from the first by cyclic permutation 
if A, LL, and 1 

There is a unique solution of this equation satisfying the boundary 


nditions in equation (2.8) given by 


P= A = 
’ J (r+a*)#(7 b*)4| 7+ s2c?/(1 +s?) |! 

dr 

B = 4 

ay | (r7+a®)i(r- b?)§[ 7 | g2p2 (1-4 s*)}! 
A 

‘ 9 9 si 9.9 - . -9\ 1a? (3.10) 
(7-+-a*)*(7 b*)![ 7 + s2¢2 (1+-s?)]? 


which may be obtained from equation (3.9) by assuming that P is a product 
fone of the cartesian coordinates and a function of A, and expressing the 
ordinate in terms of A, x, and v from equation (3.2). The constants A, B, 
may be determined from the boundary condition on the ellipsoid since 
ve have automatically satisfied the condition as x, y, and zo. The 
juations for A, B, and C' may be determined quite easily if we make use 
{the fact that if G is a function of A only, then on the ellipsoid, where A = 0, 


r oG G 
| ie (3.11) 
g* CY bh? CL 
- 9 IG CG a oG y , 0G (1 t it all (3.12) 
dr cx a cy bh? Cz sc? 
Further, for simplicity, let us write 
dr h,(A) (3.13) 
— f 3.15 
(r7-+-a*)? (7+ 5") [7 1 g2¢2 (1+ s*) |! 
nd h,(9) h,, etc. (3.14) 
Then 
rs 2(1-+-s? ys 2(1-++-s?)! 2(s*+-1 2(1-++-s?)# 
A~—{h,- -)4+ Beh, \ "\ 6 Nh, — y) 
a’ abe b2\ ° abes c2s : abes 


a on se (241779 2H) 3.15) 
Be 
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so that 
2(1-++s?)!\ | ~ 
As| h, -——_——_] + Bh, = — U(1+s?), (3.16) 
abcs ‘i 
492) - 
Ah, -8B(hy— =" “ ") reer), (3.17) 
abes 
9(] -1.92)3 = 
and o(he- = | - —sW. (3.18) 
abcs 


We see, therefore, that the motions parallel and perpendicular to the axis 
of rotation are independent, and we shall consider them separately througb- 
out the rest of this work. We may now invert equation (3.10) with respect 
to s to obtain a formal solution for P, and we may use equations (2.9)-(2.1] 
to obtain expressions for @, 6, and #@, viz.: 





SAh,(A)+ Bh (A) 2 [ Ax | By | Cz 
—T-8  —~I+s}tpe® tpt peter 5| 
.___ StA+b)+yA+@) ag 
© (Aa?)t(A+6?)HA+se2/(1+82) AR? 
sBh,(A)—. Ah (A) 2 [ Av , By , Cz oe 
ite  —1+e]\ta! Ate Adee pe 
-__ syfA+b4)—2Ata%) gy 
(A+a?)!(A-+ 6?) A+ s2e2/(1+- 8?) [HR 
and 
- _Ch;(A) Ax. me, Cz 
s jA+a?®" A+B? * A+s7e2/(1+s?) 
2sz (3.21) 


/ (A | @2) b(d re b?) {A+-s?c? (i s?) Ry’ 
where F\ is the derivative, with respect to A, of F defined in equation 
(3.3). The formal solution of the problem is now complete and it only 
remains to apply it to the problems mentioned in the introduction. 


APPENDIX A 
In this appendix we consider the branch points of A qua function of s. If, when 
$ = 8, A = Ag, then since F(s,A) = 0, near s = s, we have 


0 (s—s)F,+(A Ao) Fi + $(8 —89)*Fyy + (8 — 89 )(A—Ag) Fy, 4 $(A—Ao)* Fat... (A. 


and therefore A is a regular function of s in the neighbourhood of (s,,A 9), unless 


F, = 0. If F, = 0 when s = so, however, then near this point 
A = Ag +[—2(8—)Fi/Fyalt_ gt + (A.2) 


and has an algebraic singularity with exponent }. If F,, = 0 at s = sy as well, then 


the singularity has exponent 4, but such a case is included in F, = 0. Branch points 
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an only occur therefore when the cubic in A, equation (3.3), has two equal roots, and 


near these roots 


F, = (s—s,)*{—2F, Fy] a: (A.3) 


Further the branch point degenerates if F, 0, in which case A is a regular function 
fsnear this point, and Ff, has a simple zero. This special case, which is of importance 
in determining the surface on which the ultimate flow is not steady, may be dealt 
with quite easily if F, = 0; for, putting @ = A in equation (3.1), we see that A = Oas 

well, so that the branch points degenerate on the ellipsoid only, when 
- 0. (A.4) 
at bf (1+8?)c4 

Hence s* is real and 0 s* # 
We shall now prove that the branch points lie on the imaginary axis of s between 
the points 8 iands i. For, as A varies, equation (3.3) defines the tangential 


pencil of quadrics ASL’ — 0, (A.5) 


vhere 1, m, n, p are the tangential coordinates of the plane 


la+my+nz+p 0, (A.6) 
1+s? = 
p> l?+-m?4 — n* 0 (A.7) 
s* 
and = a*]? +-b*m?-+-c?n?—p? = 0. (A.8) 
The condition that a point (x, y,z) should lie on the quadric envelope (A.5) is clearly 


given by equation (3.3). The tangent cones through an arbitrary point Q to the 
tangential pencil of quadrics also form a pencil, and the four common planes of any 
two cones of this pencil are the tangent planes from Q to the common tangent 
levelopable of the quadric pencil (A.5). Hence, if these two cones touch, a tangent 
of the tangent developable, and therefore a generator of the tangent developable, 
passes through Q. The condition for this is that the discriminant of the cubic in A 
hould vanish, which in turn implies that the cubic should have two equal roots. 
We conclude therefore that the tangent developable of the pencil (A.5) is the surface 
n which (3.3) has two equal roots for given a, y, z, s. Further details of this geo- 
trical argument may be found in (7) in which the dual result is proved. 


Now the tangent developable of © and ’ is the envelope of 


1+-s*)4(la + my) +-2[ —s2(l2+-m?) ]# + [(a2l? + b2m?)(1 +s?) — 8?c?(1? + m?)]* = 0, (A.9) 
d hence a ty pical line of the de velopable passes through the point 

lf a*(1+-s?)— s?c’ 
l s-)2) (art b?m?)( l Ss“) ( 29%]? m?)]4 ; 


m{ b?(1-+-s*)—c?s?] : 
iiss sl te nn A.10 
s?)4{(a2l? j b?m?)(1 1 8s?) c2s?(]? Jj m*)]*’ 0) ( ] ) 


and its direc tion cosines are 


a2 1] s2 13 
i| : ie m| a : . (1+<8?)#. (A.11) 
1? +-m? 1? +-m? 
Since these lines must be real, J, m, and —s* must be real and also 0 < —s? < 1 s0 


that the branch points all lie on the imaginary axis of s in the strip —1 < ime < 1. 
Now the surface formed by eliminating J and m from a typical point on the line 
given by (A.10) and (A.11) is an equation of the eighth degree in 2, y, z, and s, and 
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of the sixth degree in a, 6, and c. In general it meets the plane z 0 in the ellipse 


a? y l 
371 1 od) 2.3 1} B21 3 = 7 (A.12 
a*(1+-s?)—c?s? © 6?(1-+-s?)—c?s? 1+<s* 
touches the ellipsoid A 1, and meets the z-axis in the four points 
l l 
2 c? a°(1 -| and 2? c?— b( 1 | ). (A.13 
\ s~ s- 
When a b the surface degenerates into the right circular cone 
(1 +-8*)(a* +y?) {isz + [a?(1 +8?) —c?s*]}}?. (A.14 


If c is very large, on the cones which pass through points near the origin we have 


so that the branch points are clustered round s = 0. 


Part IT. Untrorm Morton 
4. In this part we apply the formal solution obtained in Part I to the 
case of uniform motion, relative to the rotating axes, of an ellipsoid along 
and perpendicular to the axis of rotation. In the first problem we take 


U=V=0 and W W/s, (4.1) 
whence A B 0 and 
: ; 2(1-+-.s?)!\-1 
( —Wth, a (4.2) 
abes 


W being a constant. It is shown in Appendix B that C is regular apart 
from branch points at s?+-1 = 0. 
In the second problem we take 


U U/s, V W 0. (4.3 
Hence C0. 
IA U(h,- = a (4.4) 
abes 
IB Uh,/s, (4.5) 
2s 4 
where I = h,h,- abel p3)h (h, +h.) + abe)" (4.6) 


We also show in Appendix B that J is regular apart from branch points 
at s*+-1 0. 

If Q is a typical function in which we are interested, such as pressure ot 
velocity, then Q is already known and 


Ve ds (y > 0). (4.7) 
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The expression for Q is very complicated, so that the complete specification 
f Q involves a formidable amount of routine computing; but fortunately 
ur main interest is centred on the ultimate distribution of Q and this is 
ite easily found. We make cuts in the s-plane from the branch points 
f Q qua function of s to infinity along lines parallel to the negative real axis 
fs and replace the path of integration in (4.7) by the infinite semicircle 
res < 0 and lines closely following the cuts. The difference between the 
two integrals is equal to the residues of Q at its poles, of which there is only 
0. There are no difficulties at A 


me, in general, at s a*, for example, 


even though @, @, and @ contain A+ a? as a denominator, because A-+-a? 


un only vanish when 2 = 0, and then the term disappears from the 
functions. Since all the branch points lie on the imaginary axis of s and 
we know that A gua function of s has algebraic singularities with exponents 

at all but two of them, we may show, as in the earlier paper (2), that the 
like t-? or t-4, 


0 also tend to zero as t 


contributions to Q tend to zero ast > « The contribution 


from the special branch points at s?+- 1 > 00, but 


the demonstration is tedious and is omitted here. In general, therefore, 


the ultimate flow is determined by the residues at the pole s = 0, and in 
the next two sections we shall evaluate them. 

Special cases arise when the branch points degenerate and when one of 
them occurs at s = 0, reinforcing the pole there. The first case occurs, as 
we show in Appendix A, only on the ellipsoid, and there /, has a simple zero 
it two values of s on the imaginary axis. Hence @, @, and @, but not P, 
have simple poles at these points so that the velocity components oscillate 
finitely. The second case occurs on C, and there @, 6, and @ behave like s~?, 
while P behaves like s~! near s = 0. Hence wu, v, and w, but not P, increase 
ndefinitely with ton C. The ultimate solution is therefore not small on 
(,and in Appendix C we consider the modifications caused by the non- 
linear and the viscous terms to the flow near C. 

A third special case occurs on the axis of C when a = b, because the 
solution for the velocity has a simple pole at two values of s on the imaginary 
‘and A. 


the velocity, but not the pressure, oscillates finitely. 


xis independently of A, B, ¢ In this case, in the ultimate motion 


5. Motion parallel to the axis of rotation 
\s we have already pointed out, the ultimate motion is given, in general, 


y the residues of the various functions at the pole s 0. Now whens 0, 


and A = A, outside C, from section 3, 


C 2Wb?a/E, 


0 inside ¢ and 


(5.1) 


where E is the complete elliptic integral of the second kind of modulus 














152 K. STEWARTSON 


(a®?—6?)*/a, from Appendix B. Hence we must consider the 


regions inside 
and outside C separately. 


(a) If b?x?+-a®y? < ah? 
then oe é 
\a2b?— b222— a2y? i — 


when s is small, from equation (3.8). 
Hence, near s = 0 


P — 2(a*b?—b?2?—ay?)t 2Ez 





a ee a ee eres g2 5 3 
a*b2s ab2 ' (s*) (9.9) 
from equation (3.10), 
wD 2E : 
=—_ — onde O(s), 
C abs ° | 
u 22 91,9 9.9 9 9 
— a (a*b2- b2x*—a®y?)-*, (5.4) 
( b2s | 
D 27 | 
and ; = — (a*b?—b?x? —a2y?)-4 
( acs A | 


from equations (2.9), (2.10), and (2.11). 
In the ultimate motion, therefore. we have 


u v W (a*b? —b?x?— a2y?)-+ a 
= : a SS - ~ ——————————— {2.0 
a*y = bx ak 
w W, and P , (2b? — b2x? — a2y?)!, (5.6) 
ak 
(b) If bx? +-a2y? > ab? 
then A = Ap, defined by equation (3.5), when s = 0. Hence P. a. and? 


are bounded in the neighbourhood of s = 0 so that P. u. and v tend to zero 
as too. On the other hand, 


a 2ab2>W dr — 
sw —> _— | — rR — ass—> VU, (o./ 
EJ ri(r7+a?)'r4 b?)! 
do 
2ab?*W dr 
and so w mee | oa, E> Ow. (5.8) 
E t*(r-+-a*)#(7-+b?)! 


Xo 


This component, which is constant on any cylinder confocal with C, 





~ OW/a?+b?—2?2—y? 
~ EB ( 


1 
2,292 a as b*x?+a®y? > a*b?+ (5.9) 
0° ac" +-a*y?— ab? 
—ab?*W 
3B ty)! 


and 





as b?x?+-@?y? > oo. (5.10) 
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Inside the cylinder C' the ultimate flow consists of a uniform translation 
with velocity W along the direction of motion of the ellipsoid together with 
a swirling motion about the axis of the cylinder. The projection of the 
stream-lines on the xy-plane are ellipses similar to the projection of C on 
the same plane. The stream-lines relative to the rotating axes are there- 
fore ‘elliptic helices’ and since none of them meet C, no fluid is transferred 


Y 


from the region inside C to the region outside C. Outside C there is only 


, shearing motion parallel to the axis of rotation singular on C and tending 
to zero as b?x*+-a®y? + oc. The whole motion is independent of z and c, so 
that the ultimate flow due to a very elongated ellipsoid is the same as that 
fora disk, although it takes very much longer to achieve, since the clustering 
of the branch points near s = 0 means that the unsteady terms die out 
more slowly. On C the ultimate velocity is singular, so that the second 
order and viscous terms may no longer be neglected. A crude attempt to 
allow for them is given in Appendix C. According to the solution given 
there the velocity has a maximum near C of the order of W(a?W /zv)*, where 
vis the kinematic viscosity. 

The total pressure on the ellipsoid is in the opposite direction to the 


motion of the ellipsoid which therefore encounters a resistance 
Z | pPn dS, (5.11) 


taken over the ellipsoid, where nm is the cosine of the angle between the 


normal and the axis of C. On A 1, Pax z and hence 
, 1 «2 ee ™ 
Z=—> 3 tore as t > 00. (5.12) 


The error term when ¢ is large is O(t-? cost) and so Z is very nearly constant 
after one revolution of the fluid (¢ = 47). If the angular velocity of the 
fluid is Q, the ultimate resistance is (87/3#)WQpab?, and varies from 
7 WQpab? when b/a is small to 4*WQpa? when a = b. It is independent of c. 


6. Uniform motion at right angles to the axis of rotation 


Again, we must consider the regions outside and inside C separately. 
(a) If x*b?-+-y?a? > a*b? 


then from equation (3.8), A = Ay when s = 0. 
Further when s 0 
(a2— b?)?abcU 
(a2— b?)2+-c2(K — E)(a2E—b?K) 
(a2— b?)2ac?2(K — E)U 


and 9sB Ried Bioeth. <-cll eo 
(a*— b?)?+-c?(K — E)(a2E— 067K) 
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from equations (B.3) and (B.4). Hence, substituting in the expression for P 
we find that, as t > 0, 


2(a?—b?)abcU 


<—_—r ms ; 
(a2?— b?)?+-c?(K — E)(a27H— 67K) 
x(a*— b?) ee sid __ __ ybe(K — E) | a | 
Jr or+ary J (r+a?)i(r+ be) 
Xo x 


being the residue of P at s 0. The modulus of the complete elliptic 
integrals is (a2—b?)'/a. We may deduce u and v as ¢t > © in a similar way 
but since they are both ultimately steady it is easier to use the formulae 
... oP - . OF ? 
lim -limyv and lim— — —limu (6.4) 
tse OX ao to OY t—>a 
which may be deduced from the equations of motion. Near the cylinder 
C, in the ultimate motion, 


av bu 
x” y 
al c(a? —b?)| b?a(a? —b?) —a*®bey( K — E)|U 
an { (a?—b?)?-+-c?(K — FE) (a2? E—b?K) |(a?+-b?—2? — y?)! (ab? + ya? ah?) 
(6.5) 
so that although both w and v are singular on C, the component of velocity 
normal to C is finite, and in fact equal to 
Uc(K E)| xb*c(a2E b2.K )—ba?y(a?- b*)| . me 
9 9\9 9 ry ’ 9’ 9” ; 9 1h on ( ’ (6.6) 
| (a? —b?)?+-c?(K — E)(a? EB —b?K ) || b4x? + aty? |}! 


At large distances from the cylinder 


P (a2—b?)abe U| (a? —b?)a—be( K — E)y| (6.7) 

as : ~ —— — : 4 
3(x?+-y*)!| (a2—b?)*+-c?(K — £)(a2@H—b7K )| 

so that as 7+ -y? > 0, wand v ~ (a?+y?)-!. Further, since s® ~aP |éz 

and ¢P/éz = 0when s = 0, @is bounded at the origin, and so w — 0 outside 


the cylinder C as t-> oo. Thus the motion outside C is ultimately purely 
two-dimensional since the stream-lines then lie in planes perpendicular to 
the z-axis and the velocity components are independent of z. The motion, 
however, depends on c, the length of the semi-axis in the z-direction, in 
contradistinction to the ultimate motion when the ellipsoid moves parallel 
to the axis of rotation. When c is zero, so that the ellipsoid becomes a disk 
moving parallel to itself, there is no disturbance, and as c — 00, 

ps 2ab2U (a? —b?)y r oe dr (6.8) 

ak—b?k J r*(r7-+a?)*(7+5?)? 
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This form for P obtained by letting t + oo and then c — 00 is not the same 
1s the solution obtained by letting c + oo first and then ¢ + 00, in which case 
P + —2Ub(a+b)y ———., (6.9) 


Xo 


is may be verified by using Taylor’s theory of two-dimensional motion in 
, rotating fluid (5). Apart from the coefficients of the integrals the only 

















Fic. 1. Stream-lines for the two solutions of flow past a circular cylinder. 


Taylor's two-dimensional solution 


Limit of three-dimensional solution 


difference between the solutions is the presence of 7! in the denominator of 
equation (6.8). A sketch of the two sets of stream-lines relative to the 
cylinder is given in Fig. 1, in the particular case when a = b. The solution 
represented by equation (6.8) is not irrotational and the velocity component 
tangential to the cylinder is singular on it. Further, the normal component 
is finite there and equal to 

bal’ 


(6.10) 
(b4x?+-aty?) 


from equation (6.6). The stream-lines are, however, tangential to the body 
and the component of velocity parallel to the generators of the cylinder is 
singular at the cylinder, and zero elsewhere. It would seem therefore that 
in problems of motion in rotating fluids it is unwise to neglect edge effects 
when considering ultimate motions without first investigating them care- 
fully. 

When a = b, so that the ellipsoid becomes a spheroid, we find that, as 


> 


P— 


(6.11) 


sin - > 


2Uc(4ax | . 4a a(r? a 
7 ? 


9 


7c?+ 16a? r r2 
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where 7? = x?+-y?. Further, when c is large, 
9 y 2 2\3 
QyUT[ . a a(r?—a?)! 
p+» —2¥ | sin-17 a er : (6.12) 
7 , Y ’ 
which may be compared with the two-dimensional solution 


Ua*y 


P 9 (6.13) 
2 
(b) If xb? ya? < a*b? 
2f22 
then A 2{_ _ aa —— 02 4+ O(s4) 
\a2b?— b?a?—a?y? * 


near s = (0. After some algebra we find that, as t > 0, 
p ., S— Seite — He EOE) (6.14) 
(a?—b?)?+-c?(K — E)(a27E—b?K ) 
from which, since u and v are ultimately steady, 
- obo — =, (6.15) 


b(a?—b®)  c(@E—BK)  (a®—b2)? + e2(K— E)(@E—bK) 


as t > oo. 


Moreover, 
i c(a?—b?)U| (a2—b? )ba—ca*®y(K— E)| 6.16 
0 ee eee 
[ (a? —b?)?+-c?(K — E)(a? E—b?K) |(ab?— 62x? —ay?)! 
as t > ©. 


Inside the cylinder, therefore, the components of velocity perpendicular 
to the axis are constant, while the component parallel to the axis, although 
independent of z, has a singularity on C. The stream-lines relative to the 


ellipsoid satisfy dx dy dz 


; (6.17) 
u—l 2 w 
so that they lie in planes parallel to 
abe( K — E)+-y(a?—b?) = 0 6.18) 
2 2\ 1 
a y\} , 
and on z= cl—— =s “1 const., 6.19) 
a®_ 6? 
which are ellipsoids congruent to A 1, inscribed in C and having their 


centres on the z-axis. 


In the particular case when a = 6 these formulae simplify and we find 
that 


—_____ (4ax—cy), (6.20) 


jccglce iat. (6.21) 
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U (4%—-y)ca? : 
" U(4e— y) v—, ast—>o. (6.22) 
(7c? + 16a?)(a?—r?)8 





Further, the stream-lines lie in planes parallel to 
mcx + 4ay 0. (6.23) 
\s an example, the stream-lines which arise in the ultimate flow are shown 


in Fig. 2, when the ellipsoid becomes a sphere. It is clear that this picture 











A 
_— Ce I i ~t 
- ee i 
—— i 




















Fic. 2, Plan in xy-plane of stream-lines relative to a sphere moving 
along 2-axis. 


of the ultimate flow pattern is physically impossible to achieve. In the first 
place. about half of the stream-lines which pass through C have to make 
sharp turns before reaching it, and a real fluid would tend to break away 
irom one side of C. Secondly, the fluid crossing C has a non-zero velocity 
normal to C' while its tangential velocity is discontinuous; a real fluid could 
not do this and would probably be swept round the outside of C by neigh- 
houring particles of fluid. Both of these phenomena have been observed 
by Taylor (3). He introduced coloured liquid in front of and above the 
sphere and found that it moved straight towards C, dividing there as 
though it had struck a solid obstacle. Some of the coloured liquid went 
round one side of C remaining quite close to it, while the rest broke away 
lorming large eddies. 


Since the theoretical picture cannot be obtained in practice as far as the 
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cylinder, the theoretical solution inside C must also be different from the 
observed flow. Taylor found that the fluid inside C appeared to move 
with it, and this is not surprising in view of the difficulty fluid would have 
in crossing C’. 

The variation of the ultimate flow pattern with z and c is of interest, 
The whole motion is independent of z, and the effect of 2c, the length of the 
principal axis in the direction of the axis of rotation, is nearly the same as 
that of a scaling factor. As c > 0 the perturbed components of velocity 
tend to zero too, in agreement with expectation, but as c > 00 the ultimate 
flow does not tend to the two-dimensional solution, although, as we shall 
see below, the correct pressure balance is obtained, in agreement with the 
two-dimensional solution. Further, the two-dimensional solution is in- 
stantaneously obtained while ours is approached only very slowly. 

On the cylinder C and on the ellipsoid A = 1 the ultimate motion is 
exceptional, for on C' it increases without limit and on the ellipsoid it 
oscillates about a non-zero mean value and with periods varying from 47 
to oo. The pressure there, however, rapidly becomes steady, the unsteady 
terms varying like ¢-!(logt)-? when ¢ is large. The ultimate form for P is 
given by equation (6.14). The ultimate force on the ellipsoid is therefore 
constant and may be found by integrating equation (2.3) over A = 1. 
The coordinates of the centre of the ellipsoid relative to O’x'y’z’ are (x’, 0, 0) 
and hence if the force has components (X,Y, Z) relative to these axes 


be(a2—b?)(K — E)U 
X as 47 pabe | 9 , - < 9 ‘ A , : 9 9\9 La'|, (6.24) 
\c?( K — E)(a2 BE —b?k )+ (a2—b?)? { 
”_ EVatE—b2K eu 
Y tr pabe | (kK ; =e I : K )c a J, (6.25) 
\c2(K — F)(a2 EH —b2K )+ (a?—b?)?} 
and Z 0. 
When a b the ellipsoid degenerates into a spheroid and we have 
, 4racl | 
X 4c7rac + iy’ (6.26) 
i “P| 16a2 mc | 4 | 
; 7c? Sects 
and } 477a*cp —— es (6.27) 
: mc“ 16a" 
As c-> 0, so that as the ellipsoid tends to a cylinder P + —Uy, the 


pressure exerted on it by the fluid tends to the two-dimensional value. 
Also as a+ oo the pressure on the disk so obtained becomes singular, 
behaving like (aloga)-! but the total force exerted tends to zero. In 
Figs. 3a and 3b we give the variation of X and Y with c(ab)~! for various 


values of a/b. 
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Grace (4) obtained the force on the sphere in uniform relative motion 


by extrapolating his series solution to infinity. He found that that part 
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force which depended on U was equal to 0-83zpa3U and 
to the direction of motion of the sphere. 
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These may be compared with 0-823apa°U and 142° respectively which 
may be obtained from equations (6.26) and (6.27). 


APPENDIX B 


We prove here that 


im hel - ee te ee 
by hg eet a ee) 3 (B. 
Hi abe(1- =i | a ~ * (abc)* 
2(1+-s?)3 
and J hs - (B.2 
abes 
have no zeros in the s-plane, cut along the imaginary axis of s, from s itos 


It is convenient to write us (1+s?)/s? and consider the %-plane cut along the negative 
real axis of 4. We consider a contour D consisting of the infinite circle, both sides of 
the cut and a small circle centre ub 0, and show that the change in the argument 
of J and J round D is zero, and hence, since J and J are regular inside D, that they 
can never vanish there. To prove this we need to know the behaviour of I and J 
when ¢% is large, when & x +07, where x is real and positive, and when ¢ is small, 


(a) ob large . 


dr 2(a7E—b?Kk) - 
As ub > o£ hy “: ‘ rir b?)3(r n a2)? - “ab*(a?— b?) ’ (B.3 
0 


where K and E are the complete elliptic integrals of the first and second kinds of 


modulus (a?—b?)#/a. 


2(K—E) 


Similarly h, >——:; (B.4) 
, “ — a(a*— b?) 
9,),4 
and, since hy +hy+hs 2 : (B.5 
- abe 
oF 
J+ ——. B.6 
ab> — 
Hence, when ub is large, 
4(a°E—b°kK)(K—E é js 
I ~ ene ts —_ ) bd - + O(w-4). (B.7 
a*b?(a?— b?)* a*b*c* 
(b) ub x+t Oi: 
Let us write hy M, t iN; ete., so that 
x ce 
Mi = 1 N, r ef — 
= J (r+a®)i(r + b2)A(r e/x)3 re eae | (7+ a2)3(7 + b2)4(r7—c?/x)? 
ely 0 (B.8 


and therefore M, > 0, N, > 0 just above the cut. Similarly M, > 0, N, > 0 there. 


9 
Hence im I = M,N,+M,N, — (M,+M,) > 0 (B.9) 
> * abcy? 
and reJ = —(M,+M,) < 0 (B.10 


on the upper side of the cut. On the lower side of the cut 


imil<0O and reJ < 0. 
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& small: 


wb\F [ dr 
I lls Cast h. —~——_——_ —-_ -- ——_ 3. 
as eaies (<a) | (7+ a*ub/c?)#(7 +b%xb/c?)t(7 +1)! (5.11) 
it) 
us? a b2 : 
log | ——)+ 0), (B.12) 
( “Ww 
dr 
h (B.13 
,; | (r+ a?)?(7+b?)4(7 +c? /b)8 
U 
7 O(us log), (B.14) 
a >) y 
Qus3 . 
nd h. O(log yw). (B.15) 
a ocla b) 
Hence on this small circle 
42/5 214 6? 
I = ™ Jog(—~) + ov) (B.16) 
abe cus 
that im J has the same sign as imys and when im J = 0, rey > 0. Further reJ < 0 
erywhere on the small circle. 
Thus, as we move round D, re J is always negative so that [arg J]p 0 and J has 


zeros in D. On the part of the contour above the real axis of 4, imJ > 0 and 
when im J 0, rel 0. Similarly for the lower part of the contour, and hence J 
ver encircles the origin, which means that [arg J]p = 0 and so J never vanishes in D. 
{n attempt was made to remove the restriction that one axis of the ellipsoid should 


parallel to the axis of rotation, and to consider the translatory motion of a general 
psoid. It may be seen from this paper that there are two main problems to be 
solved. We can easily express the solution as a Bromwich integral, but before we can 


educe the ultimate flow we must show that all the branch points of the integrand 
iron the imaginary axis of s, and that it has no poles except possibly at s 0. 
first problem may be solved on the lines of Appendix A, but the second, which 
generalization of the problem of this appendix, has proved to be too difficult. 
If it is true, as seems likely, then the ultimate flow is similar to that obtained in 
his paper. 


APPENDIX C 


In all the proble ms so far considered on motion in a rotating fluid we have found 


it on C, the cylinder circumscribing the ellipsoid and having its generators parallel 
axis of rotation, the velocity components are singular in the ultimate motion, 
though this cylinder is not a solid boundary. Hence we cannot neglect viscosity 
nd the second-order terms in any discussion of the actual motion near C, and in 


s appendix we propose to make a crude estimate of their effect when the ellipsoid 
legenerates into a sphere and further consider only the singularities when the sphere 
ves parallel to the axis. When the sphere moves at right angles to the axis, the 


lestion is much more difficult not only because the velocity components are functions 


fthe azimuth angle @ as well as of r? x? +y* and z, but also because fluid is being 
transported across C. However, the exact theoretical behaviour of the fluid near C 
sacademic in this case, because a real fluid cannot turn through the sharp angle to 


ich it as required by the theor: 


In the motion « 


f a sphere of radius a parallel to the axis, if q,, gg, and q, are the 
M 


5092.29 
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radial, azimuth, and axial components of the perturbed velocities relative to ¢} 


sphere, then, when a—7 is small and positive, 


W ( 2a )’ 


Wr qz 0 and qo a a = r (C.] 
while if y—a is small and positive 
W ( 2a \3 
q; %=—90 and g,=- ~~ —} (C.2 


a \r—a 
according to the linearized theory. In the discussion of the flow near r = a we shal 
assume that q, oP /éz 0, which is certainly true in the linearized equations of 
motion and is probably nearly true in practice. Hence from the full equations of 


motion, retaining only the most significant viscous terms, we have 


cq: oq. 
q : i —- (C.3 
C2 or* 
CFa C79 
and q i v 4 (C.4 
or* 


The viscous terms are only important in a region a—8 < r < a+-8, where 8 is th: 


thickness of the ‘boundary layer’. Outside this region we may neglect them and 
then from equations (C.3) and (C.4), both g, and gg are independent of z in agreement 
with the linearized theory. If now we suppose that outside qg and q, are given by 


(C.1) and (C.2), and inside we crudely linearize equations (C.3) and (C.4), we find 


that ina oO « yp « a 0; - satisfies 
rey cq O""q- (5 
, y— (C.5 
27 \0 C2 cr" 
with boundary conditions 
’ W (2a\3 , 
q: 0 when ?r a—o and gq | 7) when r = a+0. 
a 
The solution of this equation is 
V (2a\} [1/ aW? \! z 
q: | = } erf { :; .) (r—a-+o) (C.6 
; 7\o 2 \v*9r*02? 


provided only that 6 is just sufficient to make the error function sensibly equal to 
unity when 7 = a+6 and 


a (47)" (C1 


va) 


satisfies this condition. According to this result the axial component of velocity 


5 


increases as r > ] until it reaches a maximum value of the order of magnitude of 
._(Wa?\* : 
LW | (C8 
> vz 
and afterwards decreases to zero when r ] s A similar result holds for v, but of 


course the variation occurs in the opposite direction. 
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THE RELAXATION TREATMENT OF SINGULAR 
POINTS IN POISSON’S EQUATION 


By L. C. WOODS* 
Re 1 29 November 1951; in revised form 12 June 1952] 


SUMMARY 


is harmonic or is a solution to Poisson’s equation, it may have singular points 
field or on the boundary at which it (a) has finite values, but has infinite 
ves, has logarithmic infinities, or (c) has simple discontinuities. This 
lescribes methods that can be adopted by computers using the relaxation 
jue when working in the neighbourhood of these points. Methods of obtaining 
irate derivatives and integrals in these neighbourhoods are also given. Four 
;s illustrate the methods and suggest that the accuracy is comparable with 
ybtained in similar proble1 without singularities. 


1. Introduction 
We shall suppose that the dependent variable 4 = d(x,y), of a two- 


limensional relaxation calculation is either a harmonic function, or more 


generally, satisfies Poisson’s equation, V*¢ k(x, y). If¢ has a logarithmic 


nfinity (which we will term a ‘type II’ singularity), e.g. the logarithm of the 
elocity magnitude in incompressible flow at a sharp corner or stagnation 
point (1), it is clear that the singularity cannot be ignored, but if ¢ is finite 
vith infinite derivatives (‘type I’ singularity), e.g. the velocity potential 
function at a sharp corner projecting into the fluid, it is possible to ignore 
the singularity and obtain approximate results (2). Discontinuities in ¢ 
type IIL’ singularities) could be ignored but only with serious error in the 
ution 
fthe computer tries to avoid the singularity due to a sharp corner by 
rounding it off, a very much refined mesh is required in the neighbourhood 
{the rounded-off corner, a point well illustrated in Southwell’s example X 
3). A method of avoiding the refined mesh and even allowing for a small 
mount of rounding off at the corner is given in example IIT, below. 
Singularities of type I have been considered by Motz (4) and Jeffreys (5). 
Motz used the appropriate series expansion for ¢ in the neighbourhood of 
the singular point to obtain relations between the values of ¢ at the points 
ijacent to the singularity. Jeffreys modified Motz’s method by using the 
theory of least squares to improve the values of the numerical coefficients 


New Z Scientifie Defi Corps, at present attached to the Aerodynamics 


(Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 2 (1953)] 
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appearing in these relations. Both these methods involve the use of special 
relaxation patterns with ‘awkward’ coefficients at points near the singv- 
larity. This disadvantage is overcome in the method of this paper. 
Southwell (3), Thom and Klanfer (6), and the present author (7) have 
considered singularities of type II. Southwell’s method, for a harmonic 
function ¢, involves the use of an auxiliary function % defined by 


yb = Pe log r, g? = gy? 

where the origin is at the singularity. He then puts ¢ = %-+- », and solves 
V?y = 0, y being a function without a singularity. The boundary conditions 
for » easily follow from those for 4, and the definition of 4s. It is obviously 
more convenient to be able to work directly in terms of ¢, as it takes some 
time to find the boundary conditions for 7 with any accuracy, particularly 
if the region contains more than one singularity (e.g. the two stagnation 
points usually found on an aerofoil). Furthermore there are some problems 
(7) in which the boundary conditions for ¢, and the ‘strength’ of the 
logarithmic singularity, m (see section 3), are not known exactly until the 
solution is determined. Southwell’s method would be tedious in these 
circumstances, as every alteration in the value of m would mean an altera- 
tion to the whole of the boundary condition for 7». 

Thom in unpublished work (see (7)) has used an approximate method of 
dealing with singularities of type III. 

Methods of calculating accurately integrals and derivatives in the neigh- 
bourhood of singularities, which have not, to the author’s knowledge. 
received attention before, are given in sections 5 and 8. Finally section 1] 
gives a brief account of the ‘1/r’ singularity of three-dimensional relaxation. 


2. Equations for the residuals 
The simplest difference equation approximating to Poisson’s equation 


V26 = k(x, y) is A’d,—a*k, = 0, (I 


adjacent mesh points labelled 0, 1,..., 4 in Fig. 1, a is the mesh size, and fy 
is the value of k at 0. 

If d;, are arbitrary values not necessarily satisfying (1) we can define a 
residual R, at point 0 by 


R, = A'do—a%o. (2) 


‘Relaxation’ is the process of reducing these residuals (generally in order o! 


size) to zero or near zero, by altering the values of ¢;. If we define an 
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operator A” by A’d, y ¢;—4¢,, then it can be shown (8) that a more 


o 


accurate solution is obtained by defining the residual by the equation 
Ry, = 3(4A4’¢,+A"¢go—6a?k,). (3) 


With equation (3) it is convenient to use the same ‘relaxation pattern’ 
3) as for equation (2), i.e. the simple one defined by 


om 1 6=1220, a 4 


Cho oy 
The principle involved here is met quite often in relaxation—we calculate 
residuals accurately but use simple and approximate relaxation patterns. 
This use of equation (3) (8) is similar to, but has advantages over, a method 
due to Fox (9). It involves a small amount of iteration which, however, can 
be avoided, at the cost of a more complicated relaxation pattern (10), by 


defining R, = 4A'$.+A"¢,—6ark,, (4) 
when 
if R. c R; : OR, 


— os — 20, 
Cho c dy ody 


~ 
ou 
or 
— 


Equations (3) and (4) have been introduced here as it is usually desirable 

to use these more accurate formulae in the neighbourhood of singularities 

where the independent variable is generally changing rapidly. In the treat- 

ment of singularities by the three authors quoted, only Jeffreys used an 

equation of greater accuracy than that involved in the use of (2). 
ed ds. 


on 


Gauss’s equation, — 
k dxdy = (5) 
Js . 

where s and n are respectively distances measured along and normal to a 
closed contour C’, and the surface integral is taken over the region bounded 
by C, ean often be used to check that the field lying within C is properly 


relaxed (see ex imple IIT on p. 176). 


3. The form of the dependent variable in the neighbourhood of the 
singularities 


Any solution of Poisson’s equation can be written in the form 
| 


where V2n k, 

and is harmonic. The present intention is to allow 7 to be well behaved, 
incorporating the singularity in %. Harmonic functions can be produced in 
many ways: in particular, both real and imaginary parts of re’’, where (7, @) 
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are polar coordinates in the (x, y)-plane, are harmonic functions. This is also 
true of any function of the form f(re’). It is then only necessary to choos 
the particular function which has the desired properties. 

For type I singularities, for example, the function 


os = (ret?)m, 
has finite values for positive m, but some infinite derivatives if m is not an 
integer. We can therefore take 
r’ m 
ub (asin m@+-B cos m8) : (6) 
a 
where « and £ are constants, as suitable for the treatment of type I singu- 
larities. Thus, in this case 


m Yr m 
d (asin m6+-8 cos m0) +. (7) 


The real part of log(re“’)” has a (logarithmic) singularity of type IT, while 
the imaginary part has the required form for type IIT singularities, having 
a simple discontinuity of mz along any straight line through the origin. 
Thus we can write 


d m log(“] 1m, (8) 
a 


for type IT singularities, and 
d mé-n, (9) 


for type III singularities. The introduction of the mesh size a in equations 
(7) and (8) makes the resulting forms conveniently non-dimensional. 

Of the three types of singularity, type I is the most difficult to deal with, 
as the constants « and £8 have to be estimated from the values of ¢ in the 
neighbourhood of the singularity, and also the form of the dominant term 
is more complex than in the other cases. The theory follows in the next 
section. 


4. Type I singularities 

[f a boundary on which either ¢ or its normal derivative is specified is 
discontinuous in direction at a point O, say, then ¢ will have a type | 
singularity at O. In particular if the angle between the two tangents to 
the boundary at O is z/m, then equation (7), in which the origin of r is at 0, 
gives the appropriate form for d. This case is shown in Fig. 1, where the 
square mesh is arranged so that the singularity lies on the mesh at O and 
points 1 and 9 of the mesh lie on the boundary. It is not essential that 1 
and 9 lie on the boundary—‘fictitious points’ (3) can be introduced—but 
it is essential to have 0 on the mesh. The special point 6 shown in the figure 


is atr a, 0 = z/m. 
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1) Cale ulation of « and B 


Suppose we have found an approximate solution by ignoring the singu- 


larity, i.e. by assuming a 3 = 0, then, referring to Fig. 1, if m < 4/7, 
there will be values of ¢ at points 0, 1, 2,..., 8. If the operators A’ and A” 


sre applied to equation (7), 7 can be eliminated since 


\'7 a°V2y a*k and A’n = 2a°V?n 2a*k. 
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' 
! 
\ 17 4 A ts} <i, 
ae 
23 18 12 19 ai 
Fic. | 
Thus, since for the points 0, 1,,..., 4, and 5...., 8, r/a = 0, 1 and v2 respec- 
tively, we find 
3 ; 3. 
\'b9 = @k ta ¥ sin(mnz/2)+B ¥ cos(mnz/2) (10) 
— oe 
y 0 n 0 
und 
3. 
22mA" d,,, 22imDezhyt-« > sin{m(2n+-1)/4}+B8 > cos{mm(2n+1)/4}, (11) 
n 0 n 0 
from which « and 8 ean be ealeulated. 
In the special case m } it is unnecessary to use the operator A’, as, in 
iddition to equation (10), viz. 
$,+4.+43+4,—4b) = ako +(1+v2)a+8, (12) 
we have by a similar calculation, 
bo +3+64+4,—4b) = aky+(1+~v2)a—B. 
From this ind equation (12) we find 
(] V2)ax 5 dp; d,) + dbs : ds } dy 44, aE, (13) 
and 5 3$(d, dy). (] 4) 
lft <m 2. dz does not exist, while if | < m < 1, neither ¢, nor ¢, 


exist. In these cases fictitious values of d for one or both points 4 and 8 
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must be found by extrapolation along lines not passing through O. Ifm = | 
the singularity vanishes unless either ¢ or its derivatives are discontinuous 
at 0. Simple discontinuities in ¢ are discussed in section 10, and discon. 
tinuities in the boundary derivatives could be treated similarly. If m > | 
the singularity can be ignored in most cases, as the first derivatives remain 
finite; otherwise values of $3, ¢4, and ¢, could be found by extrapolation 
for use in equations (10) and (11). 

When a and f have been computed from (10) and (11), they are used, as 
shown below, to calculate corrections to the residuals at the points neigh- 
bouring the singularity. These corrections are then relaxed, which changes 
the values of q)...., d,, then a reapplication of (10) and (11) leads to new 
values of « and £8, and so on. This iterative process (necessary only with 
type I singularities) converges rapidly, as even the initial solution, obtained 
by ignoring the singularity completely, will yield values of « and 8 reason- 
ably close to their true values (see example I in section 6). 


(6) Calculation of residuals in the neighbourhood of the singularity 


Consider point 2 in Fig. 1 for example. Applying the operator A’ to 


equation (7), we find 
P zs ; r m Z r m 
A’d, = A’ng+a0A sinmé| + BA’||—) cosmé@| , 
a 2 a 2 
where (see Fig. 1) A'd, = by) +5 +49+4¢;—4¢., 
F r\m 
, ° | . . 9 . ° 
A sin mo] 2!”(sin }ma-+-sin 3mz)+ 2” sin }mz—4 sin }mz, 
a ° ‘ 
with similar equations for 
r) m 7 
A’n, and A ( cosmé6| . 
a ‘ 


Now since A’, = ak, it follows from the above equation for A’d, that the 
residual at point 2 should be defined by 


(7) sin mo ea’|(”) cos m8) (15) 
a : a ° 


instead of by equation (2). The wand f terms in (15) are the corrections to the 
residual mentioned above. If the residuals defined by equation (2) have 


R, = A'$.—a?k,— a’ 





been relaxed, then it only remains to relax these correction terms. 
By asimilar calculation we can show that the residual defined by equation 
(3) should be replaced by 


R, = }(4A’$.+A"$.)—a2k,—3 44 | (*)" (asin m6-+-B cos mé)| — 


- x"|(7)" (asin m@+ 8 cos mo)| 
a 


9 


\ (16) 
r 


9 
« 
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m | where Ab d, +-b44+4,5+$3—4o, etc. 
rete Residuals at other points in the neighbourhood of the singularity can be 
li 10 - ~ . . 
— culated similarly. However, unless very great accuracy is required, the 
} ] . . . 
deal rections to the residuals due to a and f need be calculated only for points 
‘emair . ° > 6 . . 
main 7) 9... 8, for, if the coefficients of « and f are calculated for a point such as 
yNatior ‘ . . — 
lation 10, they will be found to be negligible. 
ed Calculation of re siduals on the boundary 
CC, as . . 
neig} When boundary gradients are specified it becomes necessary to calculate 
rem residuals on the boundaries. By a trivial extension of an equation given by 
anges : | - : 
ape he author (8) in the case V*r 0, it can be shown that if V?n = k, the 
> with ormal derivative, at point 1 for example, is given by 
ained » ON ; 9 9 10 o 2]. Olas 17 
ASO! ba . Ne Mist 457 =Not =Ng n— 3a*k, + O(a"). (17) 
n On 
Differentiating equation (7) we find 
a Ch ? : ‘ , 2On 
A’ to m| (a cos mO—B sin m@) + : 
r oo a r 00 
whence on y QO. 2 0. 
Cy Ch ay\m 1 
"pa a ~Ma ; (18) 
cy cy a 
ndin particular at point 1 we have, using (17), 
Cah . 
6a — 6ma+ not 413+ 495+ 29+ 2N9- 10y,—3a?k,. 
CY 
Eliminating » by using (7) we have finally 
ithe | 5am, |. 
“ oy }(P2+-13+ 465 + 269+ 269— 106,)—a*k, + 2ma— 
(15 "sin ma+sin 4m7+4 X 2!" sin mz)— 
8(52” cos mz+ cos 4mz+ 4 » 2!™ cos }mm 12x 2™—10), (19) 
the - 
ave Z, say. 
We can thus define a residual at point 1 by 
I : 
t10n 
, ch 
R, Z 2a . (20) 
cy 
\t point 9 and other more distant points on the boundary it will usually 
16) quite sufficient to define residuals by the form corresponding to (19) but 
gnoring the « and § terms. 
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If boundary values are specified, and it is necessary to calculate norma 
boundary gradients, then, when the relaxation is complete, equations 
similar to (19) can be used for this purpose. 

[t is useful to notice that the coefficients of « and f are symmetrical and 
antisymmetrical respectively about the line 6 = z/2m. The following table 
constructed for the particular examples calculated in section 6, gives some 
idea of the magnitude of these coefficients, when the residuals are defined 
by equations similar to equation (16). 


Point I 5 2 


2 3 
coefi. of x | as 9 ee "O01 30 "0955 0257 0955 
B) ; cs ae 0553 ° 0553 
ri. 1 ‘0431 *O1DS5 I45I O45! *2045 

BY - 1447 o185 ° 


2047 "0447 


5. Integration and differentiation near Type I singularities 
We shall first describe simple methods of integration and differentiation 
in any region of the field. To calculate { y dx, where y is any function 
known at the mesh points, a useful formula for many problems is Simpson’s 
‘one-third’ rule, i.e., referring to Fig. 1, 
: 
| x dx = ja(x,+4X0+x3)+O(0"). (21 
1 


A formula for calculating derivatives due to Bickley (10) can be written 


124(“) 4D¥%+- D¥4+- DY 2a(£ (¥*4)] L O(a>), (22 
CY} ow /0 


where Dj = ¢,—¢,, DY = $;—¢,, Dj = $1:—¢3, ete. Equation (22) is 
generally more convenient than the one-dimensional form, viz. 
| 9 Cc d 8 Du f, d O 5 (23 
aa S149 + P12 19 + O(a"), 
\OY 0 
which is, however, sometimes necessary near the boundaries. If equations 
(21) and (22) are combined we find 


36 


dx = },(18D¥+ 8D¥+ 8D¥+ Dy+ D's), (24 


me 


although in practice it is generally better to calculate the derivative an¢ 
integral separately, particularly if there are singularities in the field. The 
method of calculating derivatives and integrals in the neighbourhood 0! 
singularities now follows. 
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€ normal] Differentiating equation (7) we find 
qQuatio) ‘ ‘ 
{ r\ < | a : /Y | : 
7 sin mé cos mM6\ a4 
: n\— | (=) (= ] | | 
re 1 | 
’ Mr U . Cc7 
ng tabi ( Jeos mé ()sin m6} 2 a<" 
e , \\a a CX 
wan (25) 
» define ) | 7 . {a 
: m| ‘(2 }sin m0 +- | }cos mo x4 
iy a) l\\a \a | | 
{/y w\. bi on 
( Jeos mo — }sin m@)B| +-a- 
\\a a, ey | 


\s examples we will calculate (é¢/éx), and (é¢/éy), (see Fig. 1), on the 


sumption that there is a singularity at 0 as described in the previous 








tion. \{ point 2? equations 25) become 
Ch : CY ) 
a| m(B sin 4mz—a« cos $m) 4 a( 7] 
Ca - " 102} » 
; as. (26) 
ntiat ads - en 
—" nd a| m(8 cos 4m7-+-asin $m7)+a 
unctl \oy]. - - ey)» 
nN psol ‘ ‘ ° = 
| By direct applic ition ot equ ition (74) 
1D3+ Dp +D%o 4(5— 6) + (m1 — 3) + (Ma— is) + P2+- 9B, 
és here 
tx 2!”(sin Jma—sin 3m) + 5*”[ sin m(4a—e)—sin m($7-+ €)|—sin mz, 
itten a 2 
ty? co ma—cos $7 
ii 5 cos m(47—e)—cos m(4n- €)| —cosm7-+l, 
29) j, Pande = tan-!($). From an equation in 7 of the same form as (22) we see 
t, since » has no singularity at 0, the terms in 7 in the equation im- 
edi itely above can be rep! iced by l2a(c n OX )o 2a3(ek OX)». The term 
23) § 12a(én/éx), can then be eliminated by using the first of equations (26). This 
culation o1ves 
latl 
Cd 
12a) LDz4- De D; 
:’ (ck ‘e 
4 (l2msin sma q)B (12m cos Jmr- p)x- 2a°(' ‘ (27) 
in . ” CX} 9 
Similarly 
ve al 
1 12a| {Dy Dy+-D 
; \CY} 2 E 
00a C y 
(12m cos }ma—q’)B- 2a°(' } , (28) 
” CY} 2 
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where 


p’ = 5!"(sin m(47—e)+sin m(47+ €)|+4 x 2)" sin }mz—sin mz, 
and 
q’ = 5!" cos m(42—e)-+ cos m(47-+-€)]+4 x 2!” cos 4ma—cos mz+1, 


Normal boundary derivatives have been mentioned in the previous 
section (see equation (19)). 

For m < 1 the first derivatives of ¢ at 0 are infinite, but if necessary 
we can calculate the first derivatives of 7 at this point. From equations 
(23) (modified to apply to the x-direction) and (7) we readily find 


C ; . 
124(2) == 8D§+4,,—¢ + (8—2™){asin m7+B(cosmna—1)}; — (29) 
7) 
while from the equation corresponding to (17) for point 0, and equation 
(7) we find 
6a) — $,4.4,4+2(b, +d9)-+4b,—10b.— 30% 
ae, stot 2($, +3) +4b2—10p9—3a7ky— 
CY] o 
—af{2'"(sin }ma+sin ?mz)+4 sin $mz-+ 2 sin mz}— 
— B{2*™(cos }mz-t-cos 3m7)+4 cos 4mz7-+-2 cosm7+2}. (30) 
The examples given above should sufficiently illustrate the method of 
calculating derivatives. 
For any three collinear mesh points at which a function f of ¢ can be 
calculated (e.g. points 5, 2, 6) there is no difficulty in evaluating integrals 


6 
of the type | f(¢) dx, but for integrals terminating on the singularity 
5 


special treatment is necessary if f(¢) becomes infinite. Consider for example 
the integral of 6¢/éx from point 0 to point 10. From the first of equations 
(25) it follows that 
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10 

and it remains to evaluate | (@y/é2) dy. From equation (21) we see that this 

0 


integral can be calculated from the values of (7/02), (€n/E2),9, and (@y/€2)p. 
The value of (@/éx), follows from equations (26) and (27), (@7/éx),»9 follows 
from (22) and (25), while (@n/éx), is given by equation (29). In the special 
case when k is constant we find, after some calculation, that 
10 
8 dy — 3(17D§+12D2-4+ 8D%,4+ Dig +-bx4—$ 
- ay 36(17D3 + 12D5+8Di.+ Die+ $11 —$9)4 


CX a 


« 


0 
+f (8—2”)sin mz—4p+ 12m x 2”-1 cos 4mz—36 x 2” cos 4m fat 


+-ge{(8—2™)(cos ma— 1)—4q—12m x 2”-1sin 4ma+ 36 x 2” sin kma}B. (31) 
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9 
the expression for the integral | (é¢/@y) dx follows by a similar type of 
0 
Jculation, although in this case the equations for éy/éy are different, e.g. 
nat], »/éy), has to be caleulated from equation (19), and (@n/éy), is given by 
previous quation (30). 
The methods of this and the previous section can be easily extended to 
necessary [he more accurate finite difference formulae for calculating derivatives and 


equations |integrals. 


», Examples of Type I singularities 








J). 94 : 
ay = Rrample I: a pote ntial proble m 
equation | The first example is a potential problem with a discontinuity in the 
undary direction of 27. We shall solve the problem (referring to Fig. 2) 
¢ F fal B 
591 |6 608) 7 64419 TOI} 12 776 i 863/9 95414 1000 
yt (30 ] 1 i 1 ° °o °o 


ethod of 


57415 59q7 624) 1 683/16 765/14 858/11 95315 _|1000 
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cularity $42 |3 553|5_— S79|10 642/27 744420 491395015 _ i000 
, ] 1[-4] «=3f63J— -1 3} ro) 
xample 
juations 
500 500 500 500 730136 _844|'7__ 949]6 ico 
D E 31/89 1g ° 
Fic. 2 
| for which 6 = 1,000 on AB, d = 500 on DE, @¢/én = 0 on BCD and EA. 
hat thie § Essentially the same problem was solved by Motz (4), but he calculated his 
| olution in a field twice the size of that shown in Fig. 2. Instead of 6 = 500 
én/éx),. | on DE he put d = 0 on the extension of BA and so he had a sharp corner 
follows | at E of angle 27. Clearly 6—500 is antisymmetric about AD, and so it is 
special | only necessary to consider the region A BCDE. 


To facilitate comparison with Motz’s solution the distance of the nearest 
mesh points from ABCD was taken to be 4a, and the residuals were calcu- 
ited from equation (2), although in this particular problem the more 
iccurate formula (3) leads to results differing only by one or two units from 
those given in Fig. 2. The figures below the values of ¢ are the amounts 
to be subtracted from ¢ to obtain Motz’s solution, while the entry to the 
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right of the ¢ value at each mesh point is the amount to be subtracted from 
¢ to obtain that relaxation solution which ignores the presence of the 
singularity at E. For this problem m = 4, and so equations (13) and (14) 
are applicable. From these equations and the values of ¢ given in Fig, 2 
we have a = 0, B = 3(730—500+ 230) = 230. For this example (m =}. 
k 0) equation (15) becomes ; 


R, = A’d.+0-2758 = A’d,+63. 


This figure ‘63’, and the corresponding figures for the other mesh points 
are shown in squares in Fig. 2. The residual on the boundary was calculated 


from (see Fig. 1) 
R, 26;+4)+%9—46,— (2 >: 2?” cos }ma-+ 2”"— 4)B 
26;+¢)+¢,—446,+9, 
as this is the boundary equation corresponding to the first approximation 
(2). Of course, at first we can only estimate an approximate value of 8 
which can be improved as the solution progresses. From the relaxation 
solution, in which we assume 0, i.e. ignore the singularity, we find 
8 = 194, and so the convergence of this iterative process is very rapid. 
Integrals of the normal derivative taken along DE, EF, and GH should 
all be equal. They give in fact the difference between the stream function ¢ 
on the conductors HA and BCD. Equation (31) becomes 
10 
a OF 36(17D3+12D5+8D%,4 Dig +$11— $y) + 65 


CL 


. 


0 
180+ 65 245. 


Continuing this integral up to F we have 





ct dy = Sai. 


Cr 


. 


: 
The integrals along GH and DE are 340 and 342 units respectively. Because 
of the large contribution of the B term to these integrals (65 for EF and 116 
for DE) it is clear that as far as integration is concerned the singularity must 
be taken into account. The contour integrals are very sensitive to inaccurate 
results in the field, and in the next example much better results are obtained 
using the more accurate formula (3). 
Example II: capacity of a square section cable 

Both Gilles (11) and Jeffreys (5) have worked out this particular example, 
but only Jeffreys made any allowance for the singularity. Fig. 3 shows a 
portion of the square cable at a potential of 10,000 surrounded by a similarly 
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situated square sheath of twice the side at zero potential; JE and AD are 
ines of symmetry. The figures below the ¢ values are ‘ 1e amounts to be 
sibtracted from ¢ to give Jeffreys’s results, but as Jeffreys used only three 
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figures in his calculation, these amounts should be divided by 10 to give a 


true comparison. 


From equations (10) and (11), in which m 2 and ky = 0, we find 
Y 4.480 and B 0. Equation (16) becomes, for example, 
R, Fi 4\ ds 1 A" dbs) 1 ()-O955a. 


the last term of which, (—428), is shown in a square in Fig. 3. The corre- 
sponding figures for the two other special points are shown likewise. 
Integration of the normal derivatives has been carried out along DE, CF, 
BG, BHI, and AT, the results being 12,793, 12,795, 12,789, 12,792, and 
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12,788 units respectively, which are remarkably close. The special methods 
were required for the last two integrals which terminate on the singularity 
If the average of these integrals is multiplied by 2 x 10-*/z (see Jeffreys) 
then, by Gauss’s theorem, the resulting value is the capacity of the cable. ( 
Thus C = 1-27914x 2/7 = 0-8143 units; Jeffreys obtained C = 0-89). 
while Gilles, ignoring the singularity, found C = 0-832. Neither Jeffreys 
nor Gilles used contours terminating on or close to the singularity. An 
analytical result, which is presumably quite accurate, was given by 
Bowman (12). He found C = 0-8144. 
Example III: a torsion problem 

Southwell (3) gives a relaxation solution, based on equation (2), of the 
torsion problem for a rectangular shaft with two keyways. Fig. 4 shows 
the author’s solution, based on equation (3), of the same problem. The 
factor k has been selected so that on the mesh shown in the figure, 


ak — 400. 


GF and AF are lines of symmetry, and on the other boundaries ¢ = 0, 
Southwell had no singularity at C as he rounded off the corner so that at 
points L and M, distance }a from C, ¢ was taken to be zero, and C was 
assumed to be an interior point of the field. This was taken into account 
by refining the mesh in the neighbourhood of C to 4a, a procedure which 
more than doubled the number of mesh points. The four figure groups shown 
in Fig. 4 are the values of ¢ obtained by the method of this paper on the 
assumption that ¢ = 0 at C. We readily find for this example « = 1,859, 
B 192, and that the increments to the residuals (as defined by (3)) due 
to the a and f terms are 189, 48, 167, and 48 at points P, Q, R, and § 
respectively. The symmetry about GF doubles the increment to the 
residual at S. 

As a check on the solution, the contour integral of the normal derivative 
of ¢ was evaluated along ABCDEFA. On AB the non-zero value of k has 
to be allowed for (see (17)) when computing the derivatives, otherwise the 
method is the same as in the previous example. Simpson’s three-eighths 
rule was used over the first three mesh intervals of BL, while for the first 
interval of AB and for DE the symmetry condition was used. Values of 
the contour integral for various sections of the contour are shown in the 
figure. Because of the symmetry there is no contribution from EFA. 
Equation (31) was used over the intervals terminating on the singularity, 
although ¢, does not exist in this case—a fictitious value of zero was 
assumed. The value of the contour integral is 16,010. From equation (5) 
we see that the theoretical value of this integral is a2k x (number of mesh 
squares within the contour), ie. 40040 = 16,000. The discrepancy 
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between these results implies that the sum of the unrelaxed residuals on 
the 37 mesh points within the contour is about 10. 
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THE 


So far we cannot really compare our solution with Southwell’s since he 
lid not put ¢ 0 at C. However, if we assume that ¢ 0 at M and L, 
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nd then interpolate for the value at C in the x and y directions, we find the 
igures 323 and 404 shown in the diagram. The average of these, 363, is 
easonably close to the value given by Southwell, 354. Now this value 
extends’ over ta and so by dividing by 8 we arrive at the figure 45 as the 
ufluence of the rounding-off on the coarse mesh. This produces residuals 
tPand R which, when relaxed, result in the field of values entered below ¢. 
{these values are added to ¢ we arrive at our final solution. This procedure 
‘, of course, somewhat empirical, but it does make possible a comparison 
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between Southwell’s method and that of this paper. The values shown ty 
the right of dé when added to ¢ yield Southwell’s solution. The maximun 
difference between our solution and Southwell’s solution at P, Q, and Ris 
only about 0-2 per cent. Elsewhere the solutions diverge as we have used 
the more accurate formula (3). The main conclusion is that our method of 
first introducing a singularity, then modifying the result by interpolation, 
has saved a large amount of computation in the neighbourhood of the 
rounded-off corner. 


7. Type II singularities 

Singularities at which the function has a logarithmic infinity cannot be 
ignored as is possible with type I singularities, but on the other hand they 
require very much less algebra in their treatment. 
(a) Singularity at an interior point 

teferring to Fig. 1 we shall suppose, as before, that the singularity is at 
the point 0, i.e. dg = 0. It will be convenient to define a fictitious value 


of do, say $0: by * : 29 
0 No> (v2) 


then in any difference equation which would normally involve ¢5, we shall 
rl 


replace ¢, by the fictitious value ¢j. Equation (8), which gives the form 
of ¢ near 0 is repeated here for the reader’s convenience: 


r 
d = mlog|—)+ 7». 
a 


Now either = t 2 ~ a*ky), 


4 

a a it PY, . 2). 99) 

when, from (8) and (32), > = }( > 4¢;,—a@kp), 33 
t 

or, more accurately, 


4 8 
No 2( 4 21 ss > i % 5a*ko), 


2 


4 8 
when 0 =2(4>¢ + > $;—2m log 2—5a*ky). (34) 


2 1 o 


These equations enable us to define residuals at the point 0 in the usual 
way, e.g. if we use (33), then 


4 
Ry Zz $; _ 465 —a*ky A'dy—aP ko, etc. 
i=1 


This result is very convenient; it means that in finding ¢% we can relax al 
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wint 0 just as if no infinity exists. At the point 2, with the aid of equation 
8), we find 

A’$o = $5+46 +419 +45 —462 = A’n.+2m log 2, 
nd so we can define the residual at this point by 
R, = A’d,.—A’n.—2m log 2 = A’d,.—a?k,— 2m log 2. 
fore accurately 


R, = 3(4A'b.+A"d_) 0Pk, — 7 {44 (og _) +4'(log ‘ , 
a} 2 2 


6 | a). 
e, R, = 4(4A’¢,.+-A’d,) = a®k,—1-1924m. 
‘imilarly R; = 1(4A’$;+A"¢,)—a2k, — 0-3221m, (35) 
nd R, 1(4A'h,+ A"b,)—a*k,+0-0058m, 


ndso on. It is usually necessary, if high accuracy is required, to consider 
ints up to a distance of 3a from 0 in this way. This is no more arduous 
han type I singularities as the formulae can be calculated very quickly. 

Singularity at a boundary point 

If boundary values of ¢ are specified, the procedure is just the same as 
bove, except that we cannot use equations like (33) and (34) to calculate 

and hence ¢§. However, from equation (8) we can find values of y on 
he boundary, except at point 0, the value at which is then obtained by 
nterpolation. 45 then follows from equation (32). 

[fnormal boundary gradients are specified we proceed as follows. From 
juation (8) we find 


Ch ai C ax c > 
a m| | La—!, a? _ m z)+az. (36) 
c 


oy i oy Ox a7 


Suppose ¢d/cy is specified on the boundary y = 0, then on this boundary, 

irom (36), ed on 

a a—. (37) 

Oy coy 

>/Cy)y may not be defined, but (éy/éy), can be readily found from the 

oundary values of é¢/éy, by interpolation. Then from equation (17) applied 
) point 0, viz 


6a 7] Ns + Net 4yot 2, +2y,— 10H ,—3a*kp, 
Uf 0 


nd equations (8) and (32) we find that 


* ; , - — 2° 
dp i I , de 1, 20, 1 2d. —ba (: ") — sak, m log | ° (38) 
\~F/0 
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Using this value of d* we can now calculate the residuals at interior points 
of the field just as before. On the boundary, say at point 1, we have 


> Cd : CY a 
6a 6a ) M13 + Not 45 +2 I+ 2y— 10, —3a*h, 
CY} cy), ' ¢ 
$i3 ds 1 44; T 245 T 24,—10¢,— 3a*k, 
—m/(4 log 24 log 5), (39 


and hence 


R, (dis t dy + 4d, t 24* . 2d, 10¢6,)— ak, 2a 


od 
— 1-1924m. (40 
oy 1 

with similar equations for other boundary points. 


8. Integration and differentiation near Type II singularities 

As examples we shall calculate é¢/ey and é¢/éx at the point 2, Fig. | 
From equations (8) and (36), 

12a(2) 12a(=") 4D2 + Dt,+ Dz, (41 
\ ex 2 OX} » ‘s 

od on| : c 
— 12a} - + 12m = 4D¥+ D¥+ D¥+m(12—log 80), (42 
CY)» CY) 5 i 


and 12a 


while, on a boundary, equations of the same type as (39) are necessary. 
As an example of integration consider 


6 
r Ch 
| da. 
‘ cy 
5 
Using equation (36) we have 
6 1 6 6 
"Od - ay ,(x ° O1 ° Or 
| dx =m | J a }4 | da = Ama+ de, 
J cy ot a J oy - J oy 
5 1 5 5 


Now, from (42), 


CT 

12a !) — 4DY+ Di+ Dj—m log 80, 
CY)» ™ 

and similarly 


CTY 
12a(~") 4D¥+- D¥+- Di, — 4m(log 8+-4 log 5), 


a, 
12a !) 4D} + D¥+- DY, —4m(log 8+-4 log 5): 
cy “ = 


6 

combining these results we find 

6 

r od m 
dx 1 


dy 36( 18D T 8D¥ T SD§+ Di, 1 Dis) ~gqil8 —4log 400—log 8). (49 
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) Example IV: mapping a circle on a square 

By symmetry only one octant need be considered. We shall suppose that 
he circle in the (x, y)-plane maps into a square in the (¢, s)-plane. Clearly 
regular stars will be avoided if the calculations are made on the square 
esh in the (d,%)-plane instead of the usual (x, y)-plane. Of course the 
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joundary conditions are given in the (a, y)-plane, but for this simple example 
there is no difficulty in deducing the boundary conditions in the (¢, )-plane. 
There are several dependent variables at our disposal, but in order to illus- 
rate the logarithmic singularity we shall take log R where R is the radius 

the (a. y)-plane. This function is harmonic in both the (x, y) and (¢, #) 
lanes, having as its conjugate function @, the direction of the vector along 
vhich R is measured. The value of m is unity. 

Fig. 5 shows the (#, %s)-plane on which one quadrant of the circle is mapped. 
{Band AC are lines of symmetry, while, taking the radius of the circle to 

unity, we have log R = 0 on BC. The solution, multiplied by 10+, is 


hown ‘n the figure. log R above the line and 6 below. The value entered 


1 is (log R)*, given by equation (34), while the entries to the right of the 
esh points are the increments to the residuals due to the presence of the 
ngularity (ef. equations (35)). The value ‘6’ was found from 

*c(log R) 
dis. 
Ch 
Special methods were necessary for the two inner contours DE and FG. 


lsnoring the singularity we obtain 0-7833 and 0-3617 for 6 at E and G 
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respectively. From equation (43) we find that the effect of the singularity 
on the inner contour is to add $(0-8473) = 0-4237 to the value of the integra 
obtained by ignoring it. Similarly 0-0022 is the necessary increment to the 
other contour integral. The theoretical value of 6 on AC is clearly }p, 
(0-7854), and so the results obtained in Fig. 5 are very good. There jg q 
singularity in the mapping at the outer corner, but this does not affect the 
relaxation as both « and £ are zero at this point. 

An example of the case when m is unknown until the solution is complete 
(see remarks in the introduction) has been published elsewhere (7). The 
problem was to design an aerofoil given the velocity distribution wong the 
chord. The dependent variable was log q, where ¢ is the velocity magnitude, 
which has an infinity of strength 7/27 at a trailing (or leading) edge of 
angle 7; tr is not known until the aerofoil has been designed, but a good 
approximation can be deduced from the given velocity distribution. This 
approximation is then progressively improved as the solution is obtained, 
Each small alteration in 7 simply results in small changes in the residuals 
at the special points neighbouring the singularity. 


10. Type III singularities 
We shall take the axis of 6 in equation (9), which is appropriate to type III 
singularities, to be along the direction 0, 2 in Fig. 1, and we shall assume 
that the boundary is the line 3, 0, 1; ¢ then has a simple discontinuity of 
amount mz at point 0. From equation (9), 7, is clearly the mean value of 4 
at the discontinuity, i.e. 
a oe 
No = 3($i0+$-0) = do, say. 
Also from equation (9) it follows that 
2]. 
prot bot+bst no—4h2 = Net 5+ 0+ No— 472 = Ake, 
and so we can define a residual at point 2 by 
R, = A’$d,—a’*k,, 
where the value to be taken for 4, is the mean, ¢3. 
Similarly R; = A’¢;—a’ks, 
whereas for point 14 
Ry, = A’dy4—a7k,4+m($7-+tan-! }—4 tan! $) 
= A’b,,—@7k,4+0-0378m. 


If the more accurate definition of the residual is employed it is easily verified, 


for example, that 


R, = }(4A'¢;+A"¢;)—a*k,—}mz. 
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It follows from equation (9) that 


f ay? 4 : 4 
CO n OY C axm Cc 
a waa Get, wh © op - ees 2 (44) 


Cx 2 Ox oy ne i 


and so by exactly the same technique as for the other singularities, deriva- 
tives and integrals can be calculated in the field. Consider, for example 
(referring to Fig. 1), 13 


| & dx. 


J Y 


9 


Using equation (44) we have 


13 2a 13 13 
J P geal . ~- 
CO Xr aa CY si co7Y 
dx =m = of | l dx = 4m> log 5+ 1 dex. 
I Y Jere gg ‘ J GY 
2 0 2 2 
Now 
13 


On de sa (=) | 4(' ") : a | 
cy * ley] CY) s cy 13) 


s(18D¥ + 8D4+ 8D¥, + Dy + Dy) + 


em 18(tan—! $—4)+- 8(407—427)—tan-!$+tan-" 1}}, 
thus combining these equations we have 


13 
" Oc 


da js(18.D¥+-8D¥+- 8 DY, + Dz + D¥,)+-0-0910m. 
ct 


An important example of these singularities is ‘6’, the flow direction in 
two-dimensional incompressible fluid motion. At a stagnation point @ is 
discontinuous, while its conjugate function—the logarithm of the velocity 


magnitude—has a logarithmic infinity. 


11. A simple singularity in three-dimensional relaxation 
In three-dimensional relaxation (13) singularities of the type m/r (at 
point sources, charges, etc.) can be handled quite simply in the following way. 


[f ¢ is harmonic we define the residual at point 0 in Fig. 6 by 
64, = A’dy, say. (45) 


Now suppose we have a point charge of strength m at 0, then putting 


bu oe, (46) 
: 
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where 7 is harmonic without a singularity at 0, and r is measured from ¢ 
we have for the fictitious value of 45, 


6 6 
do No é > Ni é > Pi—m, (27) 
i=] t 
2// / 





No 


™~ 

















using (46), Now 


11 
A’, = 3 4,+-60—64, 


i 


is we 


. l : 
A’n, + m| _ 6], from (46), 


and so we can define the residual at point 1 by 


Similarly R, = A’¢,-4 2| oo Le Jm, 


v2 V3 V5 
and so on, for other neighbouring points. 


12. Conclusions 
tion can be performed in the neighbourhood of singularities, for a dependent 


field. Once the effect of the singularity on the relaxation residuals at 
neighbouring mesh points has been calculated there is no further difficulty 


of the method of this paper over previous work. 





It has been shown that relaxation, numerical differentiation. and integra- 


variable satisfying Poisson’s equation, as accurately as elsewhere in the 


the calculation proceeds in the usual way. This is an important advantage 
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It is sometimes an advantage to introduce singularities into the boundary 
conditions where none really exist. This can be done at points where rapid 
rates of change of boundary direction would otherwise demand a very fine 


mesh, e.g. at the rounded-off trailing edge of an aerofoil. 
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THE EFFECTS OF SHEAR FLEXIBILITY AND 
ROTATORY INERTIA ON THE BENDING VIBRATIONS 
OF BEAMS 


By R. W. TRAILL-NASH¥F and A. R. COLLAR 
(The University, Bristol) 
[Received 29 November 1951] 


SUMMARY 


The paper gives a fairly complete theoretical treatment of the problem of the 
lateral vibrations of a uniform beam when shear flexibility (but not shear-lag effects 
and rotatory inertia are taken into account ; frequency equations and modal forms 
are given for the nine cases defined by the various combinations of the three usual 
types of end support. Considerable attention is directed towards the higher frequen- 
cies; one matter of unusual interest is the existence of a second spectrum of fre. 
quencies, which has not previously been noted. Formulae are given by means of 
which it is possible to estimate the effects of shear flexibility and rotatory inertia as 
corrections to results obtained from studies in which these quantities are ignored. 
Two numerical examples illustrate the appearance of the second spectrum. A second 
part of the paper describes an experimental investigation conducted on a family of 
box beams; excellent agreement is obtained with frequencies calculated by th 
methods of the first part. This investigation shows the considerable importance of 
shear flexibility and the relative unimportance of rotatory inertia for a box beam 
of modern design. The paper concludes with a short theoretical study, by means 
of the ‘segmented’ approach, into the effects of flexibility in shear on the bending 
vibrations of the wing of the Bristol Brabazon I aeroplane: here also the effect is 
important. 


1. Introduction 
IN the past it has been the usual engineering practice, in the calculation of 
the natural frequencies of the lateral vibrations of beams, to neglect the 
effects of shear distortion and of rotatory inertia. The practice can be 
justified for slender beams, at least for the graver modes of vibration, since 
it is well known that in such cases the effects mentioned are small. 
Heretofore this simplified approach has been universally employed in the 
study of the vibration characteristics of aircraft and aircraft components. 
such as wings, fuselages, and propellers. However, since in recent years it 
has become necessary to discuss higher orders of vibration, and also in view 
of the changing trends of design—particularly the introduction of stressed 
skin construction in metal and the unorthodox planforms and aspect ratios 
associated with high-speed flight—an investigation into these effects was 


+ Now of the Aeronautical Research Laboratories, Fishermen’s Bend, Melbourne. 
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thought to be necessary. It was with this in mind that the present study 
was undertaken. 

The paper falls naturally into three main parts. In the first a fairly com- 
plete study is made theoretically of the case of the uniform beam, the 
appropriate frequency equations being derived for all the usual end con- 
ditions; the modified equations are also given in the case when either or 
both of the parameters, proportional respectively to shear flexibility and 
rotatory inertia, are neglected. When both parameters can be neglected the 
solutions reduce to the classical results quoted by Rayleigh (1). In this 
theoretical study the usual engineering assumptions of the Bernoulli—Euler 
bending moment—bending curvature relation, and the shear force—shear 
distortion relation, are adopted. One matter of unusual interest (academi- 
cally at least) which emerges from this study is that, at sufficiently high fre- 
quencies, a complete new spectrum of natural frequencies appears when 
both shear flexibility and rotatory inertia are taken into account; this new 
spectrum does not appear to have been noted previously. A reasonable 
physical explanation of this effect is that it is due to resonant interaction 
between the rotatory inertia forces—rotation resulting from bending motion 

and the shear stiffness. 

In the second part of the paper we discuss some experiments which were 
devised to check the calculated natural frequencies of a free-free beam for 
which the effect of shear flexibility is considerable in the second symmetric 
mode. 

In the third part a brief description is given of the results of a study of the 
effect of introduction of the shear flexibility on the bending vibrations of the 
Bristol Brabazon I aeroplane wing; the effect is very considerable. 

A possibly important limitation to the theoretical treatment of this 
problem given here is the omission of shear lag effects, which are certainly 
present in the vibration of any box beam type of structure (2). Coinciden- 
tally, the beam studied experimentally was, for other important reasons, 
so designed that shear lag was suppressed; thus the experiments conform 
with the theoretical treatment. However, it is recognized that in future 
some study of shear lag effects, in addition to those discussed here, may be 
required for a completely satisfactory solution of this problem of vibration 
of box beams. 


2. Nomenclature 

The dimensions of the quantities listed here are given in terms of mass, 
length, and time, denoted as usual by M, L, T respectively; dimensionless 
quantities are denoted by a unit. 
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Definition 

Arbitrary constants. 

Bending stiffness. 

Shear stiffness. 

Distributed moment, per unit length. 

Distributed force, per unit length. 

Arbitrary constant. 

Radius of gyration of beam section. 

A specified length. 

Applied moment. 

‘Slave’ parameters dependent on 4, «, f. 

Applied shear force. 

Xatio of angle of bending to angle of shear. 

Time. 

Length coordinate. 

Lateral displacement. 

Component of lateral displacement due to bend- 
ing. 

Circular frequency of vibration. 

B/CL*: parameter proportional to shear flexi 
bility. 

k?/ L*: parameter proportional to rotatory inertia. 

x/L: lengthwise position variable. 

y/L: lateral displacement variable. 

z/L: variable for bending component of lateral 
displacement. 

wL*,/(o/B): frequency parameter. 

Mass per unit length of beam. 

Amplitude of shearing angle. 


Part | 


BEAMS 


3. Derivation of the differential equation 


Consider an element of a uniform beam shown in Fig. 1. 


The beam is 


originally straight and lies along the axis Ox: the element is of length 6z. 


Let Y(a,t) be the total lateral displacement of the section, parallel to Oy, 


and let Z(x,t) be the component of this displacement due to bending. Let 


the shear force be Q and the moment /, reckoned positive in the direction 


Oy and the sense from x to y, respectively. If the shear stiffness of the 


section is C, the angle of shear will be Q/C; if this angle is added to that 
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arising from bending, namely, the slope ¢Z/éx, we have as the total angle 


eY @% .Q 
2M. (1) 
Ox Cx ( 


Further, if B is the bending stiffness, the Bernoulli-Euler curvature formula 
as eZ M 


ox? B° 
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Now suppose that the shear and moment arise from a distributed applied 
normal load f(x,t) and moment F(x,t), each per unit length. Equilibrium 
considerations for the element give 


a6 
od 
gM peo. (4) 
Cx 


Differentiation of (1) and (4) with respect to x, and substitution from (3) for 
Q and from (2) for Z yields 


ey M 1 : 
ex? CB — oh ”) 
¢ 2) c 7 . 
es. : ie : F — 0, (6) 
ox* Ox 


and elimination of M between (5) and (6) gives 


ay 1, 1eF_1é 
éxt BB 


to 


= ae f (7) 
B ex C ex? 

Equation (7) defines the displacement Y in terms of the applied force and 

moment distributions. When these arise from inertia loads. that is. are 

reversed mass accelerations, we have 


, o?Y 
f(z,t) = —o—_, (8) 
ot 
F(a, t) —ok? “al: . (9) 
ot*\ Cx 
From (9) we have, using (2), (5), (8) in succession, 
oF e2/M eeY o ey 
— = —ok? k? — a (10) 
ea e sale) * “ales o a 
Substitution from (8) and (10) in (7) finally yields 
oy (7 ok?\ (a o (Fe ok? a] Q (11) 
ext = \C Blades) "B\ee® " C ett} 


Equation (11) is the differential equation governing the general motion of 
the beam under inertia loads arising from its own motion. 
When the motion is sinusoidal we write 


Y (x,t) = y(a)sin wt: 


then (11) becomes an ordinary differential equation in x 
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It is convenient at this stage to write (12) in dimensionless form; we there- 





fore define €é=2/L 
n= y/L 
‘ft an atin }. (13) 
x= B/CL* 
B= 2/2 
where L is some specified length. Equation (12) then assumes the form 
ae $*(a+B) TPF a8 0. (14) 


In this equation, therefore, ¢ is a number proportional to the frequency of 
vibration; « is inversely proportional to the shear stiffness, that is, vanishes 
with the shear flexibility; and f is proportional to the rotatory inertia. 

With « zero, (14) corresponds to the equation given by Rayleigh (1); with 
,and f zero it is the equation principally studied by Rayleigh. An equation 
corresponding to (14) has been deduced and examined in relation to a simply 
supported beam by Timoshenko (3, 4); he gives a correction for the shear 
effect in the fundamental mode. Davies (5,6), following Timoshenko’s 
analysis, has also given an equation corresponding to (14); he applies it to 
the particular case of a loaded fixed-free bar, and obtains a numerical esti- 
mate for the correction due to shear on the fundamental mode of a slender 
rod loaded at the free end. An equation corresponding to (14), but with 
3 zero, has also been studied by Jacobsen (7, 8) in relation to the vibration 
problems of buildings. Goens (9), like Davies, studies the vibration of a bar 
(in this case free-free) as a means of determining Young’s modulus experi- 
mentally; he also follows Timoshenko’s analysis and notes that there will 
be a critical frequency above which the hyperbolic functions of the classical 
solution become, in this more general case, trigonometrical functions, and 
remarks that a change in the vibration mode (Schwingungstypus) may be 
expected. Kruszewski (10), like Goens, studies the simultaneous (as distinct 
from the individual) effects of rotatory inertia and shear; he limits his 
analysis to certain aspects of two of the nine cases given in the present paper 
(namely, the cantilever and the free-free beam) obtaining results equivalent 
to those of equations (39) and (49) given here. Other investigators have also 
studied particular aspects of the problem: a short bibliography of related 
papers is given following the references. 

However, in none of these investigations is the appearance of a new 
spectrum of vibrations noted; the studies are mainly restricted to the graver 
modes and to one or two particular cases of the end supports. Moreover, 
these earlier studies relate principally to solid bars. In a box beam the shear 
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stiffness is almost entirely due to the webs, and the bending stiffness maink 
to the flanges; both components, together with diaphragms, etc., contribute 
to the inertia. There are thus no simple relations between the bending sti. 
ness, shear stiffness, and rotatory inertia of the kind that exist for solid bars 
and in what follows the parameters « and f are regarded as independent, 


4. Boundary conditions for the differential equation 
Since equation (14) is of the fourth order, its general solution wili involve 
four arbitrary constants, which will be determined for particular cases by 
application of the appropriate boundary conditions. For any particular 
case, it will usually be possible to specify, at one or more points on the beam 
one or more of the following quantities: displacement, total slope, slope due 
to bending only, moment, shear force. To determine the equations for these 
quantities we proceed as follows. 
From equations (5) and (8) we have 
M @#Y ao #yY 
B as! CH 
ae , M, dy . ow? 
and if M = Mysinwt 0 yoy 
B dx? ( 
In dimensionless form the latter equation is 
LM, d*y 
B dé 
Similarly, from (4), (9), and (1), 


qd an = | ok? ail: = c) 


Cx ot*\ ox C 


1 h2ax7. (15 


or, if @ = Q)sinwt, 


oka? dM, dy 
Qo - 9 gk2y? 
7" GC dx = dx 
On substitution from (15) for 1, this may be cast into the dimensionless 
form . 
L?Q, dy dy ' 
°(1—¢?a8) —- — p*(a+ g (16 
B - dé («+B) dé 


Finally, from (1) we have, if Z zsin wt 


dz dy Qs, 


dx dx C 
and substituting for Q, from (16) we may write this in the dimensionless 
form ; ; 
dé P dy , °° dy ~ 
> (1—4?a8 xy —§+ (1+4202 ‘ (17) 
dé pe PO) 


in which € = 2/L. 
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S mai In general, it is only possible to specify any of the quantities discussed 
ntribute | above when they are constrained by some external agency to vary in a given 
ing sti. | manner, or when in constrained or unconstrained motion they vanish. It 
lid bars. | is conditions of the latter type that we shall apply in what follows. The five 
ndent types of condition may thus be summarized as: 
A. Displace ment zero: n 0. 
a ~ dr 
InvO B. Total slope zero: ] 0. 
Cases dé 
irtiey C. Slope due to be nding only (dZ dé) 2ero: 
ie be 1] 2 
d?y 2 dy 
lope di X he t (] + hx?) : = ©, 
a dé dé 
OT the SF 
d*y 
D. Moment (.,) zero: + p*a7 Q. 
ag* 
- ss d3, ° dr 
E. Shear (Q,) zero: + ¢*(a+B) te @, 
dé dé 
5. Solution of the differential equation 
The differential equation (14) is a linear equation of the fourth order, and 
thus has four solutions of the type 
uf] Ae%, (18) 
Substitution from (18) in (14) gives the characteristic equation 
At+-4?(a+ B)A?9—g?(1— a8) = 0, (19) 
which is a quadratic in A*. We may therefore factorize (19) in the form 
(A?-+- p?)(A?—q?) 0, (20) 
where 2p? = g%(a+B)+4/{b"(a—B)2+ 447}, (21) 
- 2¢? $?(a+B)+4/ {¢4(a—B)*+ 467} (22) 
— ind we note that 
{ -] Ove t < > 9 > 
p?—q? = (a+) (23) 
| peg? = $(1—$*ap). (24) 
It is to be observed, from (21), that p? is always positive; accordingly 
there are in all cases two roots A Lip which lead to components of the 
solution ; 
n = A,sin p£é+A, cos p&. 
nless We next observe, from (24), that since p? is positive, g? has the same sign 
as (1—¢?a8). Thus for those frequencies for which ¢?a is less than unity 
7 we shall have two roots A t-q, which give components of the solution 
» = A,sinhgé+ A, cosh gé. 
5092.22 
O 
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However, as the frequency of oscillation rises, ¢?a8 must ultimately exceed 
unity, and then q? will be negative; if we then write 
f= —#, (25) 
the hyperbolic functions above are replaced by 
n = A,sinr€+A,cosré. 

This last form for the components is of special interest, since (in some 
cases at least) it introduces a complete new spectrum of natural frequencies, 
as indicated in section 1. This will be exemplified later. 

For simplicity we shall adopt as the standard form of the complete 


solution n = A,sin pé+ A, cos pé+ A, sinh qé+ A, coshgé. (26) 


The modifications necessary when ¢?a8 > 1 can be obtained by writing 
q = ir in the formulae derived from (26). 
In the sequel we shall also require the following quantities and relations, 
derived from (21)—(24): 
2(p*—¢*a) = 2(9°+ 4B) = —4*(a—B)+4/{¢*(a—B)?+4¢7}, (21) 
2(¢°+ da) = 2(p?—¢?B) = ¢*(a—B)+./{b4(a—B)?+4¢7}, (28) 
and (p?—¢?x)(q?+ 47a) = (p?—¢?B)(q?+¢78) = ¢?. (29) 
Equations (27) and (28) show that the four quantities p?—d¢*x, p?—¢°8 
q’°+¢7a, q?+¢78 are always positive. 
6. The derived boundary conditions 
When the expression (26) for 7 is substituted in the five conditions A to F 
of section 4, and use is made of the various relations given in section 5 
between the ‘slave’ parameters p and q and the quantities 4, «, 8 on which 
they depend, the conditions may be written 
A. Displacement zero (n = 0) 
A, sin p£+ A, cos p€+ A, sinh gé+ A, cosh gé = 0. (30) 
B. Total slope zero (n' = dyn/dé = 0) 
p(A, cos p§— A, sin pé)+-9q(A, cosh gé+ A, sinh gé) = 0. (31) 
C. Slope due to bending zero (¢' = dt/dé = 0) 
q(p*—¢"a)(A, cos p§— A, sin pé)-4 
+ p(q?+-¢?x)(A, cosh gé+ A,sinh gé) = 0. (32) 
D. Moment zero (M, = 0) 
(p?—¢?a)(A, sin p€+ A, cos p£)— 


—(q?+ ¢?x)(A, sinh gé+ A, cosh gé) = 0. (33) 
E. Shear zero (Q, = 0) 


q(A, cos p§— A, sin p£)—p(A, cosh gé+ A,sinhgé) = 0. (34) 
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7, Simplified cases 

We shall consider, as special cases of our more general study, the simpli- 
fications which arise when 

(a) rotatory inertia can be neglected (8 = 0), 

(b) shear stiffness can be regarded as infinite (a = 0), 

c) both shear effects and rotatory inertia are neglected (a = B = 0). 
The parameters p and q in these cases have simpler forms and interrelations. 
We find from section 5 the following results. 

a) Rotatory inertia zero (8B = 0). If we write p,, q, for p, q in this case, 
we see from (27), (28), (29) 

1.2 9 °..2 9 4.2 2 
2( pi b2x) 2q7 -p?a+-4/(p4a + 4d ) | (35) 
(pi—¢?a)(qi+¢?a) = pigi = ¢° 
It is to be noted that both p? and q? are now always positive; 7j does not 
change sign as ¢ increases. 

(b) Shear stiffness infinite (x = 0). If in this case we write po, q2 for p, q, 

we have similarly 
Ip? — 2(q2-+-d2F 22.1 |(h4Q2-1 44,2 
2p = 248+ 9°B) = $B + Viet 4g") “tt 
P393 = (ps—4"*B)(G3+¢"B) = ¢ 
Again p3 and q5 are always positive. 
c) Rotatory inertia zero and shear stiffness infinite (x = B = 0). If here 


we write ps, q3 for p. 4. p= gQ=¢. (37) 


If we represent the general case, in which bending, shear, and rotational 
inertia are all taken into account, by B+-S+ RI, then the special cases (a), 
b), and (c) reduce respectively to B+-S, B+ RI, and B. The special results 
ofthis paragraph may be inserted, if necessary, into the boundary conditions 
Ato E of section 6; however, we shall here derive the frequency equations 
ind modal forms for the general case, and then state the simplified equations 


ibs) special cases 


8. Frequency equations and modal forms for various end conditions 

In the cases analysed below, except that of the cantilever, we shall regard 
the length LZ employed in the earlier analysis as half the length of the beam; 
we shall then specify the appropriate conditions at the centre (€ = 0) and 
t one end (é |). In the case of the cantilever, L is the total length, with 
conditions specified at the root (€ = 0) and tip (€ = 1); thus the cantilever 
can be regarded as one half of a beam of length 2Z clamped at the centre. 

We shall deduce, for the various end conditions given below, fre- 
queney equations and modal functions for both ¢’a8 < land ¢?a8 > 1. In 




















196 R. W. TRAILL-NASH AND A. R. COLLAR 


presenting the results we shall avoid the use of tangents, since in computa- 
tion the large numbers which may be involved are troublesome; and a sine 
will often be associated with its argument in the form, for example, p sin p 
or (sinhq)/q. A function such as (sinh q)/¢ can be expressed as a series in rk 
and there is thus continuity through the value unity to (sin r)/r when ¢ be. 
comes negative and is written as —r?; this also is desirable in computation, 

8.1. The cantilever: B4+-S+RI (¢?a8 < 1) 

At the root the displacement is zero; but since there is shear strain, only 
the component of slope due to bending vanishes. At the tip the moment 
and shear force are zero. Thus 


Tan 


- = 0, yo £° a= @, 
é=1, My = Q = 9, 
and on insertion of these conditions into the formulae of section 6, equations 
(30), (32) to (34) give 
A,+A, = 0, q(p?—¢a)A,+ p(q?+¢?a)A, = 0, 
(p?—¢*x)(A, sin p+ A, cos p)—(q?+¢?x)(Azsinhq+A,coshq) = 0, 
q(A, cos p— A, sin p)— p(A, coshq+ A,sinhq) = 0. 
The condition that not all the quantities A should vanish is that the deter- 
minant of their coefficients vanishes. Expanding the determinant and using 
the relations of section 5, we obtain 
sin p sinh gq 
Pog 
Any value of ¢ which satisfies this equation accordingly gives a natural 


2+-{2+¢?(«a—B)*}cos p cosh g—¢?(a+ 8) = 0. (38) 


frequency of vibration. 
The corresponding displacement function, if K is an arbitrary constant, is 
ten ( PS pla? +-d2a)sin pe —a(p*—d*a)sinh gE} + 
p q 
+{(q?+¢?x)cos p+ (p?—¢?x)cosh q}(cos pé—coshqé), (39 
in which it is to be remembered that ¢ must be a solution of the frequency 
equation (38). 
For B+S+ RI (¢?a08 > 1). When ¢?a8 > 1, gq? becomes negative; if we 
write g? = —r?, (38) becomes 
sin p sinr 


P r 


2+{2+¢4?(a—B)*}cos p cos r—¢?(a+ B) 


and the modal function is 


0 (40) 


Ky es Race “|e —x)sin pé—r(p?—¢?x)sin ré} 
p r 
Lf 


+{(7r?—¢?x)cos p—(p?—¢?x)cos r}(cos p§—cos r€). (41) 
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For B+ S. We put 8 = 0, and using the results of section 7 obtain the 
frequency equation 
2-+-(2+-47a°)cos p, cosh q,—d¢asin p, sinhg, = 0 (42) 
and the modal function 
= > hh (pj sin p, —qj sinh q, €)+ 
P3 41 
(pj cos p,+4qj cosh q,)(cos p, €—coshq,&). (43) 
For B+RI. We put « = 0 in (38) and (39); the frequency equation 
— 2+ (2+¢°8")cos p, cosh qg—¢B sin p, sinh q, = 0 (44) 
ind the modal function 
Ky = (q28in po— pe sinh qy)(q2 sin py €—p, sinh q, €)-+4 
(q3 cos p.+ p cosh q,)(cos py &—coshq,€). (45) 
For B. With « = B = 0, the frequency equation and modal function 
wssume the familiar forms 
1+cos vd cosh vd = 0, (46) 
Ky = (sin vé—sinh v¢)(sin €véd—sinh év¢)-4 


(cos \é+ cosh vd)(cos Evd—cosh Ev). (47) 


8.2. Free-free beam, symmetric modes 

In this and in the other types of support considered below the results are 
rather simpler than those of the cantilever, since two of the constants A, in 
each case vanish. 

With the origin (€ = 0) at the centre of the beam and with € = 1 at one 


end, the boundary conditions are evidently 


§ 1) 7’ Vo 0. 
i, M Qo 0. 


The first two conditions yield, on substitution in (31) and (34), 
pA, T qAy, 0 | 
qA,—pAy 0) 
which, since p?+-q? + 0, require A, = A, = 0. The remaining two con- 
ditions are then given by (33) and (34) as 
(p? b*x )A, cos p (q? + ¢?x)A, cosh q QO, 
gA,sinp+pA,sinh¢g = 0. 
From these two equations we now obtain the following results (in all cases 


the frequency equation is given first and the modal function second). 
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For B+S+RI (¢?aB < 1) 
p(p?—¢?a)cos p sinh q+-9(q?+- ¢?a)sin p coshg = 0, (48) 
Kn = qsinp cosh gé—psinhq cos pé. (49) 


For B+S+RI (a8 > 1) 


p(p?—¢?a)cos p sin r—r(7?—¢2x)sin p cosr = 0, (50) 

Kn = rsinpcosr€—psinr cos pé. (51) 

For B+S q, COs p, sinh g,+>/p, sin p, coshq, 0, (52) 
Ky = q,8in p, coshq, €—p, sinh q, cos p, €. (53) 

For B+ RI p.cos p, sinh q,+-43 sin p, cosh q, = 0, (54) 
Ky = q28in p, cosh gq, €—p, sinh q, cos py €. (55) 

For B cos vd sinh vd-+ sin vd cosh vd = 0, (56) 
Ky = sin vd cosh €v¢—sinh v¢ cos Evd. (57) 


8.3. Free-free beam, antisymmetric modes 
Here the centre of the beam is a nodal point, and since there is symmetry, 
the moment must vanish. We therefore have 


é = (), n V 0. 


a 
é= 1, My Qo = 9. 
The first two conditions show that A, and A, vanish: the remaining two 
conditions give the frequency equation and modal function: 

For B 7 S+ RI (¢?a8 <x 1) 
p(p?—¢*«)sin p cosh q—q(q?+ ¢2a)cos psinhg = 0, (58) 
Ky = peoshqsin pé+ ¢ cos p sinh gé. (59) 

For B+S+ RI (¢?a8 > 1) 
p(p?—¢?x)sin p cos r—r(r?—¢?a)cos psinr = 0, (60) 
Kn = peosrsin pé—r cos psin ré. (61) 
The reduction to the simpler cases in which shear or rotatory inertia or 
both are suppressed has been exemplified earlier and need not be given here. 


8.4. Simply supported—free beam 

Since the conditions at the centre of the beam in the case dealt with 
in section 8.3 are precisely those of simple support, the results give the 
frequencies and modes of a beam of length LZ simply supported at one end 
and free at the other. 
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8.5. Simply supported beam: symmetric modes 

The end conditions for this case are 

& 0, y = @, 0, 

g L, n = My = 0. 


The first two conditions give A, = A, = 0; the frequency equation derived 
from the remaining two is simply 


cos pcoshq = 0. 
Thus we have 


COS p 0, (62 


— 


Ky = cos pé, (63) 
and for B+S RI (d¢?a8 1), 
cos p cosr 0. (64) 
This is satisfied either by the vanishing of cos p or cosr; the latter solution 
clearly gives a new spectrum of frequencies. 
The modal functions are 
Ky cos p& (cos p 0) | 
Ky cos ré (cosr = 0) J 


As before, the derivation of the simpler cases need not be given. 


8.6. Simply supported beam: antisymmetric modes 


Here the conditions at the central nodal point are equivalent to those of 


a simple support; hence 
& 0. ” M, 0. 


€é= 1, 9= A, = @, 


and the frequency equations and modal functions are as follows. 
For B+S RI (h7a8 : 1), 


sinp = 0, (66) 
Kn = sin pé. (67) 
For B+S+ RI (d2a8 1), 
sinpsinr = 0, (68) 
Ky sin p& (sin p = 0) \ (69) 
Ky sin r€ (sinr = 0) J 


It may be remarked here that since we have a condition of simple support 
in the centre, this section gives in alternation both the symmetrical and 
antisymmetrical modes and frequencies for a beam of length L. 
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8.7. Beam encastré at each end: symmetric modes 
The boundary conditions are 
P= 0; n' V5 == (). 
Ee = 1, n Cc’ = 0, 
and we obtain the following frequency equations and modal functions: 
For B+S+RI (d?a8 < 1) 
q(p?—¢?«)sin p cosh q+ p(q?+¢2a)cos psinhg = 0. (70) 
Ky = cos p cosh gé—cosh q cos p€. (71) 
For B+-S+ RI (d? xB > #) 
r(p*+-¢?x)sin p cos r— p(r?— ¢a)cos psinr = 0, (72) 


K n COs p cos r€ — COS r COS pé. 73 


8.8. Beam encastré at each end: antisymmetric modes 


The boundary conditions are 
E= 0, n M, = 0, 
4 - ” C 
and for B+ S+RI (2x8 < 1) 
p(q’+ $°a)sin p cosh q—q(p>—¢?a)cos pginhg = 0, (74) 
Ky = sinhq sin pé—sin p sinh gé. (75) 


i 
and for B4+-S+ RI (d?a8 > 1) 


p(r?—¢?a)sin p cos r—r(p?—d2a)cos psinr = 0. (76) 
Ky = sinrsin p§—sin psinr€. (77) 
8.9. Beam encastré at one end and simply supported at the other 


The conditions here (if Z is the length of the beam) are precisely those 
of section 8.8, which accordingly gives all the results for this case. 


9. Nature of the motion in the new spectrum 

In section 8 are given all the frequency equations and modal functions 
for any combination of the three usual types of end support, namely, rigidly 
built-in, simply supported, and free. Clearly other types, such as a yielding 
support in which the displacement is in a given proportion to the shear force, 
could be applied without much difficulty. The equations are given both for 
¢*a8 < land > 1; and in certain cases it is immediately apparent that an 
additional new spectrum of frequencies appears when ¢2x8 > 1. 

We shall consider briefly some aspects of the nature of the motion in this 
new spectrum. It has been suggested earlier that the new spectrum may be 
in the nature of resonance due to interaction between the shearing forces 
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nd the rotatory inertia forces, with the rotation arising from bending. We 
navy therefore profit from an examination of the relative contributions of 
ending and shear to the motion in the two spectra. 

The current angle of shear, from (1), is y’—’, the angular displacement 
lue to bending being ¢’. Denoting the ratio of angles of bending to shear 


vy R. we find from (17) 





E vn” + (1+ 2a2)n’ 
C xn” —P? ala +-B)n’ 


0 v(1 t){n” +47(a+B)n’}+(1 *aB)n’. (78) 


The interpretation of this equation, except in special cases, is a matter 
ff some complexity. In the case of simply-supported beams, however, the 
two spectra separate completely (see sections 8.5, 8.6) and we find for the 

irst spectrum _m 97 

n P q > 

'so that, by (78), R is independent of spanwise position, and is defined by 
0) b2a(1 t){h?(ax+B)— p*}+ 46?(1—¢? a8). 

Using (23) to (25) this reduces to 


> 


$2a(1+ R) = p?. (79) 


Now the natural frequencies in this first spectrum correspond to given 
values of p (defined by sin p 0 or cos p 0); we therefore require to 
eliminate d. For this purpose we use equation (19), writing — p? for A*. and 
we obtain the quadratic equation in ¢?a8: 

(d2xB)? (d2aB)51 + p*(a t B)} t piap Q. (80) 
Clearly the same equation holds if we write r for p; in fact (see (21), (22)) the 
lower root of (80) gives the value of ¢ corresponding to a given value of p 
lorexample, p = 47); the upper root gives the value of ¢ corresponding to 
the same value of r. Elimination of ¢2 between (79) and (80) now yields 
p?a(1+ R)?—f1+ p%(a+p)\(1+ R)+ pp = 0. (81) 
This equation also has two roots; from (79) it is apparent that for a given p 
or r) the larger value of (1+ R) corresponds to the smaller value of ¢?, that 
is, to p; and similarly the smaller value of (1+ R) corresponds to r. We may 
therefore write the solutions of (81) in the forms 


R {l—p?(a—B)}+.,/[f1 p?(a— B)}?+- 42a] 


4r2a| 
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These two equations show that in the first spectrum, since R,, is always 
positive, the bending and shearing angles are in phase; on the other hand. 
R, is always negative, and hence in the second spectrum the motions are in 
antiphase. Again, f is usually small compared with a, and hence we find 


approximately R, = (p*a)-1, eg 


‘ 
It follows that in the first spectrum the proportion of bending in a given 
motion decreases steadily as p increases, and the contribution of shear 
increases appropriately. In the second spectrum the contributions of bend- 
ing and shear are roughly equal but in opposition; this implies that there 
may be strong rotations of the sections and strong shearing action while the 
lateral displacement is slight. 

The case 7 = 0 requires special consideration. When 7 becomes vanish- 
ingly small, R, tends to the limit —1, which implies that the total lateral 
displacement vanishes. The same result might be inferred from (69), which 
is formally satisfied by r = 0. However, these results have been deduced 
from a differential equation in which the lateral displacement is the depen- 
dent variable, and a solution giving zero lateral displacement everywhere 
requires closer examination. 

If in (14) we put a8 = 1 (corresponding to r = 0) the equation reduces 
d*y 1 Sa | B d?n ss 
dé" of dé ‘ 
and for general values of «, 8 the only solution of this equation for simply 


supported end conditions is 7 = 0 at all points. However, this does not 
necessarily imply that no motion is possible. If6is the amplitude of shearing 


to 


angle, and @ is used as the dependent variable in place of 7 in the analysis 
of section 3, the same differential equation results; but in this case the 
equation, with ¢2a8 = 1, admits a solution @ = constant (+ 0). The exis- 
tence of this solution is demonstrated physically by use of equations (1) to 


(4) with Y = 0, and hence f = 0; the equations are then 
0) - inl A (82) 
Ox Cc 
OL . M (83) 
Cx? B 
e@ _o. (84) 
Ox 
G4: M F — 0. (85) 
ae 


From (84) we observe that Q is a constant with respect to x; hence from 
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99) eZ/éx is constant, and from (83) M = 0. Thus (85) becomes, on use 


f (Q) 
i \ } 


ae 0 Sait () = 0 
OX OU"\ CX 
nd with @Z/éa @ sin wt this yields 
(C—ok?w?)é = 0, 
that @ is not zero provided C = ck®w?; (86) 
hat is ¢7a8 = 1 
Thus the motion corresponding to r = O is a pure shearing oscillation of the 


am, with the rotational inertia forces at any section in equilibrium with 
the shearing actions, and with zero lateral displacement. 

It may be noted that the results just established (for example, (86)) are 
ndependent of L, and accordingly may relate either to a simply supported 
yam or to an infinite beam. 

The existence of a solution at r = 0 involving pure shearing motion, and 
the earlier considerations concerning relative phases, lead to the suggestion 
f section 1 that the new spectrum is due to a modified resonance between 
rotational inertia forces and the shearing stiffness forces in the beam. 


10. Relation between the vibrations of an infinite beam and the 
higher modes of finite beams 


The vibration characteristics of finite beams depend, of course, on the end 
onditions, as illustrated in section 8. However, the influence of the end 
conditions decreases as the order of the mode increases—that is, as the wave- 
length decreases—and the modal form of the beam approaches that of a 
portion of an infinite beam vibrating with the same wave-length. Thus a 
study of the infinite beam serves to give an indication of the effects of shear 
nd rotatory inertia on the higher modes of a beam with any end conditions. 

To begin such a study we note that the fundamental mode of anti- 
symmetrical vibration of a simply supported beam of length 2Z is identical 
with that of a corresponding portion of an infinite beam, since the conditions 
it the nodal points are those of simple support. Now the fundamental 
solution of equation (68) for either p or r is p 7 or r =z. Inserting 


this value into equation (80) gives 


2467 ah }] T 7° ( x B)! T V[{f1- 7*(a-+-B)}*—474ap]. (87) 


The positive sign for the radical (yielding ¢2a8 > 1) gives ¢,, while the 
hegative sign gives ¢,,; since for r = 7 and p = a, ¢, > ¢,. When either « 
or f 


p becomes vanishingly small, the lower root survives, while the upper root 
tends to infinity. 
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In assessing the effects of shear and rotatory inertia in the higher modes, 
therefore, we must restrict ourselves to a comparison of the lower root of 
(87) and the solution derived from it when « and B vanish. If we make y 
vanishingly small, the lower root of (87) yields 

d? a/(1-+ 7B); (88) 


finally, when f also vanishes, we may write 


Accordingly, in general, 


(¢-) _{1+77(a+8)}—,/[{1+-27(a+B)}?— 42408] 
op ‘ 

The ratio ¢/¢, is plotted against « in Fig. 2 for a range of values of 8. 
Some care is necessary in the interpretation of this figure. Suppose, for 
example, that we wish to estimate ¢/¢, for the sixth natural frequency of a 
free-free beam (both symmetrical and antisymmetrical modes). There will 


(89) 


be seven nodal points,+ and a complete wave-length will occupy the space 
between three: that is, about one-third of the length of the beam. 

Thus if a, B/CL? and B, = k?/L* have been evaluated for the beam, 
based on L, the half-length, the values of « and B to be used in (89) must be 
about 9a, and 98), since they relate to a wave-length of about one-third of 
the beam length; and in general, for an assessment of ¢/¢,, in the nth natural 
frequency, the basic values ay and 8, must be multiplied by about (4n)* to 
obtain the values of « and £ to be used in Fig. 2 


11. Some considerations relating to rotatory inertia 

In practice the effect of rotatory inertia is frequently very small compared 
with that of shear flexibility. Even for the solid rods discussed by Timo- 
shenko (3,4), the shear effect is four times the rotatory inertia effect, and 
with box beams the rotatory inertia can become quite insignificant. It may 
therefore often be omitted without serious error, and a very considerable 





simplification of the problem is thereby achieved. However, if it is desired 
to estimate the rotatory inertia effect as a correction to a study in which it 
is omitted, we may use the formula, derived from (88) with and without 8 
dpert_ Ree (90) 

PB y (1+7*B) 

This formula is derived on the assumption that « = 0. The justification 
for this step lies in the fact that Fig. 2 shows the effect of variation of 8 
to be greatest when a = 0. This result may be expected on the physical 
grounds that, in a vibration mode of given amplitude, the component of 


+ We are here considering only the first spectrum, since dz corresponds to « = B = 0, 
when the second spectrum is absent. 
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displacement due to shear increases with a, and in consequence the compo- 
nent due to bending becomes less. The rotation of the sections thus de. 
creases as « increases, and the rotatory inertia effects will be smaller. Thug 
(90) should give an upper bound to the correction for rotatory inertia. An 
example of its use is given in Part IT. 

The ratio ($p..p7)/¢g is given by the top curve of Fig. 2 if 8 is regarded as 
the abscissa. The remarks of section 10 relating to an effective value of 8 
for a harmonic also apply here. 


12. A numerical example 

In Part II some experiments with a built-up box beam are described; for 
this beam B has the value 0-00325, while « can be varied between about 0-08 
and 0-20. To provide a simple but not unreal numerical illustration, « and 
B have been selected to conform roughly with these values: we take 

x+B = 0-200, 
x8 - 0-0006; 
so that « = 0-19695, B = 0-00305. 

In Fig. 3 are given curves of p, q, r against ¢, for 0 < ¢ < 100, and for 
the above values of « and 8. The change from q tor occurs for 0-0006¢? = 1, 
ord = 40-8 nearly. It will be seen that the curve of p is virtually a straight 
line passing close to the origin; it is sensibly coincident with its asymptote 
p = dv\a = 0-4444. Similarly, r rapidly approaches the asymptote 

r = dvB = 0-058¢. 

In the same figure are also given curves of p,, q, obtained using the same 
value of « but with 8 = 0; p, is indistinguishable from p, while q,, after 
early coincidence with q, rapidly tends to the asymptotic value 1/va = 2-252. 

These relations have been applied to determine the frequencies and modes 
for two of the cases of section 8, namely, the symmetric modes of the 
free-free beam and the antisymmetric modes of the simply supported beam. 
The results and conclusions are summarized briefly below. 


12.1. Free-free beam (symmetric modes) 
The natural frequencies for this case are determined by equations (48), 
(50), and (52). Solutions have been obtained by regarding the left-hand 
sides of the equations as functions of ¢ and plotting against ¢ to determine 
the zeros. Table I shows the values of ¢ corresponding to natural frequencies. 
It will be observed that, in the range covered, there are two natural fre- 
quencies—corresponding to values of ¢ of 51-1 and 94-8—appearing in the 
complete solution which are absent when rotatory inertia is neglected: the 
comparison is facilitated by the fact that, apart from the appearance of 
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these frequencies, there is little difference between the two cases. These two 
frequencies belong to the second spectrum. Although they do not separate 
out, as in the case of simple support, they may be identified in several ways, 
In the first place, they interrupt the steady (almost arithmetical) progres. 
sion of the natural frequencies; secondly, the modal functions are unusual: 
thirdly, the number of nodal points (given in the table) also interrupts the 
usual progressive increase. 


TABLE I 


Values of ¢ for Natural Frequencies 








B+S+RI1 4°2 | 12°5 | 20°r | 27°4 | 34°6 | 41°4 | 47°7 | 51°1 
B+sS 4°2 12°5 | 20°1 27°4 | 34°7 | 4I°9 | 49° 

Number of nodes I 2 3 4 5 6 ri 7 
(half-beam) 

B+S+RI1 56°7 | 63°5 | 70°5 | 77° | 84°7 | 91°5 | 94°83 | 99°2 
B+S 56°3 | 63-4 | 70°5 | 77°6 | 84:7 | 918 98°9 
Number of nodes 8 9 10 II 12 13 13 I4 


(halt-beam) 

The displacement functions » for the natural frequencies corresponding 
to ¢ = 47-7 and 51-1 are shown in Fig. 4. Both have seven nodes in the 
half-beam; the first has the type of displacement at € = 1 usually associated 
with a free end, but the curve for the new spectrum shows an unusual con- 
figuration at the free end. 

Some other features in Fig. 4 deserve remark. At the free end (€ = 1) 
there is, of course, no applied moment or shear force. Equations (15) and 
16) therefore give 48 
a) 8 0=> 7 fxn, 

0 = 7”"+¢7?(a+B)y’. 


In the simplified theory, which assumes « and f to vanish, these conditions 
would require that the curvature »” and its derivative 7” should be zero; 
in this more general case, however, neither vanishes. 


Again, equation (78), on use of (51), can be cast into the form 
K(n'—f') = ¢a(sinr sin pé—sin p sin ré), 
with Ky’ = p*sinrsin pé—?r* sin psin ré. 


These two functions vary differently across the span, so that it is not pos- 
sible, except in an arbitrary manner, to effect any comparison between the 
contributions of bending and shearing actions with the total slope in the 
two cases: this difficulty was mentioned in section 9. 
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12.2. Simply supported beam (antisymmetric modes) 

For this case we employ equations (66) and (68). The second spectrum 
now separates out completely; the frequencies and modal forms are wp. 
related to those of the first spectrum. Since p and p, are sensibly the same. 
we may construct a table analogous to Table I as follows (Table IT): the 
first row gives the values of ¢ for which sin p or sinr vanishes, while the 
second gives the zeros of sin 7. 

TABLE I] 


Values of ¢ for Natural Frequencies 











> | 
B+S+RI 58 | 13°4 | 20°7 ‘9 | 35°0 





27 4o'8 42°1 19°3 
B+S 5S | 13°4 | 20°7.] 27°99 | 35° 42°1 19°3 
Number of nodes I 3 5 7 9 (o) 11 13 

(excluding ends) 
B+S+RI 56°4 | 63°5 | 702 | 70°6 | 77:7 | 848 | 91-9 | 9Q9°0 
B+S 56-4 62:5 70°6 he 84°5 gI°9 99°00 
Number of nodes 15 $7 I 19 21 23 25 7 
(excluding ends) 


Once again the frequencies in both cases are roughly in arithmetical pro- 
40-8 (r = 0) and 70-2 (r =z). As 


)) gives no lateral displacement, so 


gression, except for the two cases ¢ 
has been remarked above, the case r 
that there are no nodal points in the usual sense. The curve of lateral dis- 
placement is identical for 6 = 5-8 and ¢ 
the sections are. of course. very different. Equation (81) in this case has 
the two roots 


R,, = 0-520, R, -0-990; 
and since R = ¢’/(y’—Z’), the rotation ¢’ is given by 
R 
2 0-342 (db = 5:8) 
n R+1 
99-0 (¢ = 70-2) 
Part I] 


EXPERIMENTS ON THE VIBRATIONS OF A QUASI-UNIFORM BEAM 
13. General considerations 

In this part we describe some experiments devised to check the more 
important conclusions of the analysis given in Part I. During a preliminary 
study of the problems involved in such experiments the following considera- 
tions became apparent. 

(a) Any attempt to check all the results by experiment would involve a 
prohibitive amount of time and labour; the scope of the experiments must 
accordingly be severely limited. 


70-2; the relative rotations of 





(b) 
tain | 
First 


ment 
that 
end ¢ 
ing ] 
sistit 
the | 
dime 

(¢, 
pair 
pairs 
thel 
dent 
stiff. 
para 

(d 
abo 
(¢). 
(In 
to 6 


mea 








cal pr 
7). As 
nent, s 
ral dis 
tions of 


] 


ase nas 


BEAM 
e mort 
minar 


sidera- 


volve 


Ss must 














SHEAR FLEXIBILITY AND ROTATORY INERTIA 211 


(b) While experiments with solid beams of rectangular section have cer- 
tain obvious attractions, they are ruled out by either of two considerations. 
First, to achieve values of « or 8 which would be significant from the experi- 
mental viewpoint, the length to depth ratio of such beams would be so small 
that the specimen would resemble an ingot rather than a beam; undesirable 
end effects might then arise. Again, such beams are rarely used in engineer- 
ng practice, being very much heavier, for a given duty, than a beam con- 
sisting of flanges and a web or webs; and in aircraft practice particularly, 
the flanges and webs may be quite thin in comparison with the overall 
limensions of the beam section. 

c) In view of (b) it was decided to construct a box beam consisting of a 
pair of flanges, a pair of webs, and a set of diaphragms; alternative sets of 
pairs of flanges and webs of different thicknesses would be provided so that 
the bending and shear stiffnesses could, to a first approximation, be indepen- 
dently varied. In this manner an appreciable range of variation of the shear 
stiffness parameter « would be provided. The plane of vibration is, of course, 
parallel to the webs. 

d) For a solid beam of rectangular section, « and f are related; B has 

bout one-third the value of «. For the practical type of beam described in 

however, 8 is unrelated to « and is by comparison quite insignificant. 
In the set of beams eventually chosen the value of « varied from about 25 
to 60 times the value of 8.) It was therefore decided not to attempt to 
measure the effects of variation of B. 

¢) Because of the insignificantly small value of 8 in practical cases, the 
second spectrum discussed in Part I appears only as overtones of high order; 
this is illustrated in section 12. Such overtones are very difficult to excite 
ind identify; it was therefore decided not to attempt this. 

f) Any set of beams constructed must be sufficiently large to render 
practicable the ready interchange of the sets of flanges and of webs without 
the use of special bolts or tools. At the same time, the specimens must not 


t 


e SO large that examination in a laboratory would become difficult. These 
considerations fixed the dimensions of the specimens within close limits. 
However, it appeared that, because of their lightness and stiffness, beams 
i these proportions would have frequencies, even in the gravest modes of 
vibration, so high that special techniques would have to be devised for their 
excitation and measurement, unless they could be artificially brought into 
ower ranges 

It was therefore decided to use massive transverse diaphragms closely 
spaced throughout the length of the beam. These diaphragms would not 
contribute t 


the bending or shear stiffnesses (except perhaps very in- 
directly) and would provide a mass loading which would very materially 
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reduce the frequencies of vibration. At the same time, the rigidity and close 
spacing of the diaphragms would imply that buckling of the webs and 
flanges, and local panel vibration or ‘panting’, would be prevented; such 
phenomena would provide undesirable complications. 

While the above considerations were the determining factors in the design 
of the specimens, there is the additional feature, associated with the rigidity 
of the diaphragms, that they virtually exclude any important shear-lag 
effects. In the theory developed in Part I the assumption is made that a 
plane cross-section of the undistorted beam remains plane during the vibra- 
tions; thus experiments on beams with very stiff diaphragms will tend to 
conform with this theory. 

(g) The inclusion of discrete diaphragms, and even the decision to employ 
bolts to provide interchangeability, implies some ‘‘eparture from the uni- 
formity assumed in Part I. However, if the number of diaphragms is large, 
the beam specimens can be regarded as quasi-uniform, and the departures 
from the results for a uniform beam are slight. This point is dealt with more 
fully in Part IIT. 

(h) In practice it is very difficult to provide a rigidly built-in or simply 
supported condition. However stiff and massive the body to which attach- 
ment is made, there is always a degree of local elasticity, the effect of which 
is difficult to estimate. Accordingly, the free-free condition only was 
selected for test. 

(i) For the range of frequency to be examined, electromagnetic excitation 
of the natural modes is more suitable than mechanical excitation. In addi- 
tion, it was decided to apply the exciting force at the mid-point of the span, 
so that only symmetrical modes of oscillation would result. 


14. Details of the beam specimens 

The beam section and layout are as shown in Fig. 5. Each beam is 36-5 in. 
long, and is made up by bolting a pair of flanges and a pair of webs to 
nineteen mild steel diaphragms, each 0-5 in. thick, 5-15in. broad, and 3-75in. 
deep, and spaced at 2-in. intervals throughout the length. In addition, the 
edges of the webs and flanges are bolted together as shown, so that the 
overall breadth is 5-75 in. 

The webs and flanges are of Alclad DTD 610; attachment to the dia- 
phragms is effected by means of 4 BA bolts in tapped holes, while 6 BA bolts 
and nuts at 4-in. intervals join the webs to the flanges. To ensure accurate 
interchangeability and freedom from structural damping all holes in the 
plates are accurately drilled and are close fits on the bolts. 

+ The writers acknowledge with gratitude the assistance of the Bristol Aeroplane Company 


in making up the specimens for test and in carrying out material tests for stiffness calcula- 
tions (see section 15). 
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Fig. 5. 

Three sets of flanges and three sets of webs were made up, giving a possible 
nine specimens having different values of the stiffness parameter. The 
nominal thicknesses are as follows: 

Webs: 0-017 in. 0-024in. 0-038 in. 
Flanges: 0-048 in. 0:°064in. 0-078 in. 
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The specimens are numbered as in Table IIT: 


TABLE ITI 


Flange 0048 in. o-064in. 0°078 in. 
Web Specimen number 
0°038 in I 4 7 
0'024 in 2 5 8 
0°OI7 In 3 6 9 


It was found in practice impossible to work with Specimen 9; the combina- 
tion of thick flange and thin web produced some buckling and ultimately 
cracking of the webs. 

[t is to be noted that, in any computations on these specimens, they were 
regarded as of effective length 38 in., each diaphragm lying at the mid-point 
of a section 2 in. long. The flanges and webs are thus imagined to extend 
0-75 in. beyond the outer faces of the end diaphragms. 


15. Determination of stiffnesses and stiffness parameters 

Values of the stiffnesses of the eight specimens studied were found both 
by calculation and test. The mean actual thicknesses of the flange and web 
plates were carefully measured, and the values of Young’s modulus for the 
various plates were also found. From these quantities the stiffnesses are 
found by a simple calculation. The shear stiffness C is assumed to be en- 
tirely due to the webs; on the other hand, while the principal contribution 
to the bending stiffness B is from the flanges, there is an appreciable contri- 
bution from the webs. No allowance can readily be made for the contribu- 
tions to the stiffnesses from the bolts or diaphragms, and the calculated 
stiffnesses do not contain any such correction. 

In view of this last uncertainty, experimental values of B and C were 
obtained by loading tests of the completed beams, as assembled for vibration 
tests, in a 10-ton Denison machine. In the main series of tests a four-point 
loading system was applied, so that the central portions of the beams were 
subjected to pure bending, while the outer portions were in addition strained 
in shear. In further check tests a three-point loading system was applied. 
From the measured deflexions at various points the stiffnesses B and C 
could be readily calculated. Since both B and C are determined from differ- 
ences of deflexions. small errors in measurement can produce relatively large 
errors in the values of B and C; in particular, an over-estimate of B leads 
to a rather larger under-estimate of C, and vice versa. It is concluded that 


the measured values of B are probably accurate to within 3 or 4 per cent.:; 


the values of C to within about 10 per cent. It may be noted that, for 


calculation of the fundamental frequency of vibration at least, the overall 
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securacy is much better than would be suggested by the accuracy of B and 
( individually: in both vibration test and stiffness test B and C are com- 
bined in much the same way, so that the errors involved in disentangling 
the individual values of B and C from the stiffness test tend to disappear 
in the recombination for the frequency calculation (see section 17). 

The measured values of B and C are compared with the calculated values 
in Table [V; on the whole, the agreement is very good, though it appears 
that there is an appreciable contribution to B from some source not con- 
sidered in the calculated values—presumably this is an effect resulting from 
the presence of the diaphragms. Table IV also contains values of « (based 
on the stiffness test measurements and on a half-length L 19 in.) and of 
the density parameter o (in lb./in. divided by g in in./sec.”) obtained by 


weighing the completed beams. 
TABLE IV 


From slifjness lests 


B eg x o 
I lb. in.? lb. units mass/in. 
°) 2 10°| o-94 10°| 0083 3°88 x 1073 

25°4 o'72 | 0-098 3°86 

2671 0°40 0-180 3°83 

35°4 0°98 0085 3°Q1 

35°0 0°04 ors! 3°90 

33°0 0°46 0°202 3°88 

$1°O I°I3 o'Iol 3°96 


0°62 0°177 3°94 


16. The rotatory inertia parameter 
In the calculation of the rotatory inertia parameter f it is assumed that 
the mass of the beam is uniformly distributed throughout the equivalent 


lencth 27, and over a rectangular cross-section of the dimensions of the 


transverse diaphragms. The radius of gyration / is thus given in terms of 
the depth d 3°75 in.) by the formula 
k? = jd’, 
and hence. with L 19 in 
3 k?/ L? 0:00325 (91) 


lor all the beams. While the different flanges and webs will in practice 
modify this value of 8 a little, it is in fact so small, and its influence so slight, 


that refinement is quite unnecessary. 


17. Vibration tests: Comparison of theory and experiment 
A free-free condition implies that there shall be no constraints on the 
beam under examination: in practice, however, the weight of the beam must 
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be supported. Supports at the nodal points of the beam are obviously the 
ideal; but, while the positions of these points may be calculated, they are 
not known exactly before the experiment is performed, and in any case, vary 
with the mode under investigation. It is therefore more practicable to use 
a suspension which is very ‘soft’ in the direction of displacement in vibra- 
tion. Such a support was achieved by hanging the beams, with the webs in 
horizontal planes, from an overhead support by four vertical wires each 
about 8-5 ft. long, attached at points along the centre lines of the flanges, 
equidistant from the centre. These wires permitted small movements in the 
horizontal plane with virtually no constraint, but prevented vertical move- 
ment or rolling about any horizontal axis. 

Vibration excitation was provided by a De Havilland electromagnetic 
moving coil unit; the moving coil was attached through the central point 
of one flange to the middle diaphragm. The vibrator unit was energized 
from either of two alternators, the frequency of excitation being deter- 
mined by the alternator speed. 

An amplitude pick-up system incorporating photo-electric cells was em- 
ployed to recognize resonant frequencies; with development, measurements 
of modal displacements could also have been made. For each beam, the 
fundamental and first overtone frequencies for symmetric vibration were 
recorded; the results are contained in Table V. 


TABLE V 


Frequencies (cycles per second) 





Fundamental ( f;) First overtone ( f») 





Calculated Calculated 
Spec.| Expt. | B+S+RI| B+S| B | Expt. | B+S+RI| B+S) B 


I 183 186 188 211 641 636 654 | II4!I 
2 170 173 176 200 550 | 572 |} 585 1084 
3 162 162 163 204 480 469 476 1103 
4 202 207 209 235 676 704 | 720 1272 
5 188 192 193 234 559 574 584 1268 
6 178 177 178 224 503 502 509 232 
7 216 219 221 254 702 716 732 1373 
5 200 198 200 249 576 576 534 1344 











Comparison of the figures in Table V is facilitated by the supplementary 
Table VI giving the ratios of calculated to observed frequencies. 

To calculate the frequencies we employ the formulae of section 8.2 in 
conjunction with the data in Table IV. Equation (48) determines the values 
of ¢ corresponding to natural frequencies; for purposes of calculation, it was 
convenient to obtain curves of variation of ¢, and ¢, (the subscript denotes 
the order of the mode) with a, for 8B = 0-00325. Corresponding curves were 
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TABLE VI 


Ratio of Calculated to Observed Frequencies 





Fundamenta First overtone 
S+kl B+S b B+S+RI1 B+S | B 
1°02 °c IIs 0°99 1°02 1°78 
I*Id 1°04 1-00 1°97 
fore) Io! 1°25 0°95 0'99 2°30 
1°02 1°03 1°16 1°04 1°06 1°38 
1° I 1°24 1°03 1°04 2°27 
o'9 1°C 1°20 100 1°Ol 2°45 
7 1-O! I‘o 1d 1°02 1°04 1°96 
9 me 1°24 100 1‘o! 2°33 
Mea O1 I 1°21 1°Ol 1°03 2°12 


lso found from (52) with 8 = 0; the isolated values on the latter curves at 
0 are, of course, the first two zeros of (56). These curves are given in 
Fig. 6; their general characteristics are very similar to those of Fig. 2, which 
relate to the infinite beam or to the condition of simple support. 
From these curves, values of ¢ are read off for the values of « obtained 
for the beams and listed in Table IV: the frequencies are then determined 
see equation (13)) by the formula 


f=asn (le (92) 


Table V gives calculated values of f, and f, for the cases B+ S+ RI, B+S 
that is, 8B = 0), and B (that is, « = B = 0). The values of f, and f, deter- 
mined from the vibration experiments are also given. 


18. Comments on the comparison 

The comparison between observed and calculated frequency contained in 
Table VI shows that there is a very good correlation between the observed 
frequencies and those calculated for the case B+ S-+- RI; for the funda- 
mental, the differences range from —1 to +2 per cent., and for the first 
overtone, from —2 to +4 per cent., the mean error in both cases being 

| per cent. We may conclude that, at least if measured values of the 
stiffnesses are employed, the formulae of section 8 may be used with 
onfidence for frequency calculations. 

Suppression of the rotatory inertia parameter 8 makes very little differ- 
ence; the calculated frequencies are all raised a little, so that the mean error 
nthe fundamental rises to + 2 per cent. and in the first overtone to +3 per 
cent. Unless 8 happens in a particular case to be considerable, therefore, 
tis probably sufficiently accurate in practice to make calculations neglect- 
ng the effect of rotatory inertia; an estimate of the importance of the error 
in readily be made by the use of equation (90). In the present instance 
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If the beams are assumed rigid in shear, however, the errors introduced 
ire very large, and (as would be expected) increase with the value of a 
neglected. An average error of over 20 per cent. is noted for the funda- 
mental; in the first overtone, for Specimen 6 (a = 0-202) the frequency 
\btained for a B 0 is nearly 2} times the correct value as obtained 
oth by calculation and observation. It is to be concluded that shear effects 
should always be included in calculations of frequency of vibration. The 
results of Part III add point to this conclusion. 

The experimental values of the natural frequencies have been reduced to 
the dimensionless forms ¢, and ¢, by means of equation (92) and are plotted 
n Fig. 6 for comparison with the theoretical curves for a uniform free-free 


eam. 


Part II] 


[HE EFFECT OF SHEAR STIFFNESS ON THE BENDING VIBRATIONS 
OF THE BRISTOL BRABAZON I WING 
19. A preliminary study 

The theory of Part I relates to a uniform beam, and the experiments of 
Part II. devised to check the theory, therefore also relate to a quasi-uniform 
beam. In considering the practical problem of the vibration of an aircraft 
wing. however, we are faced with the treatment of a non-uniform structure; 
it is usually impossible to establish—let alone solve—a differential equation 
for such a structure, since the spanwise distribution of mass and of bending 
ind shear stiffness have no analytical form, and frequently contain discon- 
tinuities. A different approach is therefore required. 

The usual procedure—which is basically what we adopt in the sequel—is 
to regard the mass as concentrated at a number of points along the span, 
these points having massless elastic interconnexions with the appropriate 
stiffnesses. This approximate representation of the true system has a finite 
number of degrees of freedom, which may be specified in any convenient 
manner; the equations of motion are then set up and solved. Matrix algebra 
is useful for the setting up of the equations, and—for example—matrix 
iteration processes may be employed for the solution of the problem (11). 

In this part this method has been applied to demonstrate the importance 
of the inclusion of flexibility in shear in calculations of the bending frequen- 

ies of the wing of the Bristol Brabazon I aeroplane; in view of the results 
of Part I, it is assumed that rotatory inertia effects may be neglected. 
\ preliminary study was, however, considered desirable; it has a double 
purpose. 


In the first instance, it was desired to obtain a check on the validity of 
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this ‘lumped mass’ or ‘segmented’ representation of the continuous struc. 
ture, and to examine the degree of approximation involved. For this 
purpose, it was decided to compare the results of calculations on the 
symmetrical free-free vibrations of a uniform beam (given in section 8.2) 
with results obtained by the ‘lumped mass’ approach. 

In the second place, it will be observed that the experiments described 
in Part II were, in fact, conducted on a physical system approximating 
closely to the ‘lumped mass’ representation of a uniform beam. The pre- 
liminary study can therefore give an estimate of the order of the difference, 
arising from this cause, to be expected between the calculated and experi- 
mental results discussed in Part II. 

Accordingly, a segmented representation of a free-free uniform beam, 
with nineteen concentrated masses as in the experiments, was examined by 
the matrix iteration process to obtain the lowest two frequencies of sym- 
metric vibration: this, with 8B = 0, is in fact a problem in nine degrees of 
freedom. Calculations were made of the fundamental and first overtone 
frequencies for « = 0 and « = 0-16 (with B = 0 in both cases) since this 
type of comparison was to be made in the case of the aircraft wing. These 
fout results are themselves compared with the corresponding results for the 
continuous beam, which are among the calculations on which the curves of 
Fig. 6 are based: the comparison is given in Table VIT. 


TaBLeE VII 
Symmetric Modes of Free-Free Beam (B = 0) 


Shear parameter x ° o'16 


4. | ¢ | ty | ds 


Values of $,, $» for 

(a) Continuous beam 593 30°23 | 4 
615 30°54 | 4° 
Ratio (b)/(a) 1'004 roro | 1 


nun 


(b) Segmented beam 





It will be seen that the error implicit in the segmented approach is in 
fact very small. The comparison is in general accord with results given 
recently by Duncan (12), who has shown that provided the masses are 
concentrated at the centres of the segments, the error is ultimately propor- 
tional to 1/n?, where n is the number of segments. We may conclude that 
the segmented approximation may be adopted with reasonable confidence 
to discuss the characteristics of anon-uniform structure such as an aeroplane 
wing. 

It may be remarked that the segmented structure has a slightly higher 
frequency than the continuous structure it represents. Thus, for an accurate 
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comparison with the experimental results (which relate to effectively seg- 
mented structures) the calculated frequencies of Table V should be slightly 
increased, to the extent indicated by Table VII. In the main, this slightly 
increases the differences between theory and experiment; but the changes 
are very small, and the comparison remains very satisfactory. 


20. Results for the wing of the Brabazon I 

The object of the present investigation was, as has been previously 
indicated, to assess the effect of shear flexibility on the bending vibrations 
of an aeroplane wing; with this in mind, it is expedient to idealize the wing 
somewhat, so that torsional vibrations are excluded. 

Curves giving the variation along the span of the mass, bending stiffness, 
and shear stiffness of the Brabazon I wing were kindly supplied to the writers 
by the Bristol Aeroplane Company: these data were applied to an idealized 
representation in the form of a straight cantilever with the mass concen- 
trated on the flexural axis at eight points along the span. The first five 
natural frequencies of this system were then calculated using the matrix 
iteration process, first for the shear stiffness as given, and then with the shear 
stiffness assumed infinite. The results are compared in Table VIII. 


TaBLeE VIII 
Bending Frequencies of Bristol Brabazon I Wing 
(cantilever) 


tency (cycles| sec.) 


{/} 
te Infinite 
Iness | shear stiffness Ratio 
I 1'092 1°165 1°07 
975 247 1°10 
10°93 1°43 
} 5 19°35 I°OI 


29°21 2°04 


It will be seen that the errors due to ignoring flexibility in shear can be 
very considerable. In flutter calculations it is frequently necessary to in- 
clude at least two modes of wing bending, and occasionally three modes: 
clearly, therefore, there is a probability of serious error if shear flexibility 


is not taken into account. 
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SUMMARY 
Extremum principles are applied to the problem of estimating loads sufficient to 
ause pronounced plastic yielding in notched bars. The theory is two-dimensional 
ind an ideal plastic-rigid material is assumed. Upper bounds to the constraint factors 
n pure bending are obtained for deep notches with circular and wedge-shaped roots, 
nd for shallow notches of any shape. The estimates for deep notches are calculated 


ins of slip-line fields constructed to satisfy all conditions in the plastic region ; 


but since the associated stress distributions in the rigid zones are not examined, it 


snot known theoretically whether these are complete solutions. A simple method is 
presented for ascertaining whether the rate of plastic work is positive in a given 
p-line field. 

Experiments are described in which wide bars of copper, stainless steel, and mild 
steel, deeply notched on one side only with a single wedge-shaped notch, were bent 
by pure couples applied at each end. The yield-point couples, the surface deforma- 

m, and the regions of plastic deformation (revealed by etching bent mild steel 

cimens) agreed closely with the theory. 
Introduction 


THE papel! 


s concerned with the estimation of applied loads necessary to 
cause plastic bending of a notched bar. Provided that work-hardening is 
negligible in the range of strain under consideration, the yield-point can be 
located as a well-defined bend on the load-distortion curve. An unnotched 
specimen, whose cross-section is the same as the minimum section of the 
notched bar, will yield at a lower value of the load under the same system 
of loading ; the ratio of the actual yield-point load to this lower load is 
defined as the constraint factor. The investigation is directed particularly 
to an estimation of constraint factors for notches of different shape. The 
problem is treated as two-dimensional, the strain being zero in the direction 
parallel to the length of the notch. This condition is approximately satisfied, 
with the exception of narrow zones near the edges, if the bar is sufficiently 
wide (at least ten times the thickness at the minimum section, as a rough 
estimate). The specimen is therefore taken to be a wide rectangular block 
notched through the thickness either on one side only, or symmetrically on 
opposite sides. A notch is described as shallow or deep according to 
whether the zone of plastic deformation does or does not spread to the 
notched side of the bar. The yield-point load of a deeply notched bar is 


[Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 2 (1953)] 
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independent of the precise depth of the notches. The system of loading 
considered in detail is bending by pure couples applied at the ends; the bar 
is assumed so long that the yield-point load is independent of the precise 
method of applying the couple. It must be emphasized that we are con- 
sidering only the yield-point in slow bending, and neither impact tests nor 
the problem of fracture are discussed in this paper. 


General considerations 

It is convenient to introduce the assumption of an incompressible plastic. 
rigid material which is rigid up to a well-defined yield-point and does not 
harden. This enables the yield-point load to be calculated with precision, 
and the result should be a good approximation for a real metal which has 
a sharp yield-point and a small subsequent rate of hardening. A detailed 
justification for this approach has been given by Hill (1, ch. ix). Briefly, 
the assumption of incompressibility in plane strain implies that the state 
of stress in any element is a pure shear combined with hydrostatic pressure, 
while the maximum shear stress is a constant k throughout the plastic 
region. This leads to an accurate estimate of the state of stress in the plane 
of flow even if the problem is not statically determined;t in the plastic 
region the stress state is defined by a slip-line field obeying the Hencky 
relations.{ The accompanying error in the stress normal to the plane of 
flow may be quite large, but a knowledge of this stress is of little interest. 
The further assumption of rigidity, which is equivalent to allowing the shear 
modulus G to increase without limit, introduces little error if the boundary 
conditions do not specify the absolute values of any displacements. The 
stress distribution under given surface loads is then independent of ¢ to 
the approximation involved in neglecting changes in the shape of the body; 
these changes are negligible below the yield-point. The velocity relations 
along the slip lines reduce to Geiringer’s equations when G is infinite. The 
problems discussed here are not statically determined even though all the 
boundary conditions relate to stresses only, and it will be necessary, there- 
fore, to consider both stress and velocity solutions simultaneously. 

The internal state of stress in a partly plastic body at the yield-point 
depends, in part, on the previous loading programme, and to obtain a 
complete solution it would be essential to follow the development of the 
plastic zone. However, Hill (2) has shown that the part of the plastic zone 
where deformation occurs at the yield-point depends only on the current 
surface tractions and not on the previous loading programme. It is sufficient, 


+ A problem is said to be statically determined when the stress solution can be calculated 


uniquely from the stress equations and boundary conditions alone, without reference to the 
associated velocity solution ; for discussion see ref. (1), p. 242. 
+ 


For an account of the theory of the slip-line field, see ref. (1), Ch. VI. 
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therefore, when only the yield-point load is required, to adopt the usual 
procedure of finding a slip-line field whose associated stress and velocity 
distributions satisfy the boundary conditions at the yield-point only and do 
not involve negative plastic work. Such a slip-line field will give the unique 
yield-point load provided that the associated stress distribution in the rigid 
zone does not violate the yield criterion. An incomplete solution is one in 
which certain of the above conditions are not satisfied or are left unexamined. 
Such incomplete solutions may sometimes be used to calculate approxima- 
tions from above or below to the yield-point load by means of certain 
extremum principles for a plastic-rigid body (2). These will now be stated 
ina form applicable to plane strain problems such as the present one, which 
involve either a single load or proportional loading, and in which there are 
no unknown constraints. 

Suppose that the body actually yields under stresses F, = Fr, over a part 
‘, of its surface S, r; being the given load ratios. Consider any equilibrium 
state that does not violate the yield condition at any point where deforma- 
tion actually occurs, and for which the corresponding boundary stresses 
we F* = F*r, on Sp. Then, by the maximum work principle for a finite 
ho (3; 2 
body (3 and 2), F > F*. (1) 
ind hence F is unique 

For the complementary minimum principle consider any velocity distri- 
bution u* having zero divergence. The velocity may be tangentially dis- 
continuous across certain internal surfaces S7,. Then it can be proved (2) 


that. for a uniformly hardened body in plane strain, 
_ ff * c . , ‘ wk < 
F [ nuk dSp < kf y* dV + | [w*] ay], (2) 


where y* is the maximum shear strain-rate (twice the tensor component) 
and |u*| is the absolute value of the velocity discontinuity across the 
surface S7, 

[f, in plane strain, a slip-line field can be constructed to satisfy all stress 
and velocity boundary conditions at the yield-point, the problem is still 
not completely solved unless it can be shown that (i) the rate of plastic 
work is positive throughout the field, and (ii) the yield criterion is not 
violated in the rigid zones. A method of calculating the rate of plastic work 
in a slip-line field is given in the following section. There is at present, 
however, no general method for deciding whether the second condition is 
satisfied. If the yield condition is in fact exceeded in certain parts of the 
rigid zone, such a solution is nevertheless valid for a material which is 


sufficiently harder in these parts. It follows immediately from (1) that the 
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surface stresses F* F*r, corresponding to this solution will tend to be an 


F : F*. (3 


over-estimate, i.e. 


Conversely, if two such slip-line fields can be constructed (one of which 
may or may not be the correct solution) so that F < Ff < F%, the yield 
condition is certainly exceeded in the rigid zone of the second field in some 


part which lies within the deforming zone of the first. 


Plastic work in the slip-line field 

Suppose that a slip-line field is constructed to satisfy the stress and 
velocity boundary conditions and we wish to discover quickly and easil; 
whether the rate of plastic work is positive at all points; more generally 
we may wish to ascertain what velocity distributions are compatible with 
a certain form of slip-line field. 


A 





\ 
os 
\YO 


Fic. 1. A curvilinear element in the slip-line field. 


Following the usual notation, the two sets of slip lines are labelled a- and 
B-lines, the «and sets being distinguished by the directions of the maximum 
shear stress k acting on a curvilinear element bounded by slip lines (Fig. 1). 
Let ¢ be the anti-clockwise angular rotation of an a-line from some fixed 
direction. It is convenient to introduce the quantities (7, 7) to define a 
point in the field. They are the coordinates of the point P under considera- 
tion, referred to axes passing through some fixed origin O and parallel to 
the respective slip directions at P (Fig. 1). The geometry of the field may 
also be described by the radii of curvature R and S of the «- and f-lines. 


Denoting by ds,, dsg the elemental are lengths measured in the positive 
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lirections along the slip lines, these radii are defined by the equations 
| l a 
R cs Ss O88 
Hence, at a point (.¢, 7) the radii are 
r cy " 
R S Y ~ ZX. (4) 
@ Ch 


Consider a 
work 1S ky dI 


lenote time-derivatives. 


where y 


small volume element dV at a point P. The rate of plastic 


is the maximum shear strain-rate at P and dots 


Let wu and v be the velocity components at P 


referred to the a- and f-directions. Expressed in terms of these velocities, 


» 


cs , R ' Ss" 


(5) 


Now it has been shown that a velocity field can be transformed into a slip- 


ine field (Green, 5). This is evident from a comparison between the velocity 


und (%, 7)-relations along the slip lines, viz.: 





du—vdd—=0 along an a-line 

dvtudd—=90 alonga f-line 
und dj+%dé = 0 along an a-line ) 

di—jdé—0 alonga f-line J 


The transformation equations which will be used here are 
u i’. v £ , db d’, x-> v’. 
Substituting for w and v in equation (5), and noting that the angular varia- 
tions in the original slip-line field and in the transformed-velocity field are 
identical, we deduce that "Re 
y - ; (6) 
S R 
where R’ and S’ are the radii of curvature in the transformed field, defined 
ina similar manner to R and S. 
If, by examining the general forms of the known (or postulated) slip-line 
R’/ R) in equation 
6) are found to have the same sign, the sign of y may be ascertained without 


and transformed-velocity fields, the terms (S’/S) and ( 


laborious calculation of the radii. For example, it can be shown by this 
method that the rate of work is everywhere positive in the solutions given 
by the writer for indentation between two flat smooth dies (5). If the 
terms are of opposite sign over all or part of the slip-line field, the shape 
of the transformed-velocity field may sometimes be altered to make the 
signs everywhere the same by superposing on the original velocity distribu- 
tion a uniform rotation about any fixed point. For suppose that the body 


is rotated in the clockwise direction with uniform angular velocity w about 
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O, the arbitrarily chosen origin of the (%, 7)-coordinates of the slip-line field, 
Then, denoting new quantities by the suffix 1, the new velocity com. 
ponents at P relative to the instantaneous slip-line directions are 


Uy U-t+wi, vy V—we. 
Hence, by (4), the new radii of the transformed-velocity field are given by 
R; R’+wR, S} S’tw/S. 


The rate of plastic work, of course, remains unaffected. If the relative 
rotation between elements in different parts of the deforming zone js 
sufficiently small, the terms (S/S) and (— R{/R) can be arranged to be of 
the same sign at all points by a suitable choice of w. The choice of 0 has 
no effect on R{ and 8}, and so a uniform velocity of translation changes the 
origin but not the shape of the transformed- 
velocity field. ; 

Velocity discontinuities are not covered by 
the above analysis and must be considered 
separately. The sense of a velocity discontinuity 
across a given slip line must be chosen so that 
the work of instantaneous shear as an element 
crosses the slip line is positive. 

Let us examine the possibility of a velocity 
discontinuity across one or more of the slip 





lines in the configuration shown in Fig. 2. 
There are two distinct pairs of slip lines through the point O, the pair of 
x- (and of f-) directions being adjacent but not coincident. O is either a 
stress singularity, or lies on a line of stress discontinuity which bisects 
the angle AOA’. Consider velocities at O, and suppose, for convenience, 
that in the region bounded by OA and OB the material is at rest. In the 
regions AOA’ and BOB’ the velocities at O must be directed along OA and 
OB for the plastic work to be positive. Such velocities, however, imply 
velocity discontinuities across OA’ and/or OB’ which involve negative 
plastic work. Hence, the velocities at O must be the same in all regions and 
no discontinuity can take place across any slip line ending on such a stress 
singularity or line of stress discontinuity. A rigid region bounded by slip 
lines through a point such as O can only rotate relative to O. 

Finally, the angle between an a- and a f-direction at a point must be at 
least $7 to avoid exceeding the maximum shear stress k on an intermediate 
plane through the point. Hence a velocity discontinuity is only possible if 
the «- (and B-) directions at O coincide, the slip lines are continuous, and 
there is no stress discontinuity. 
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Pure bending of a notched bar 


Suppose that equal and opposite pure couples, M, are applied to the ends 
fa bar of unit width notched either on one side only, or symmetrically 
n both sides, the thickness of the minimum section at the notch being a. 
For notches of all shapes and any depth an upper bound to M is obtainedt 
from a deformation consisting of the rigid rotation of each end by shearing 
slong two equal circular ares of angle 2A, say, spanning the minimum section 
ff the bar to form a rigid pivot; (2) gives an upper bound of ka?A/2 sin*A to 
MV. and the smallest value of this is 0-69ka?, when tan A 2A. An obvious 
ower bound is provided by the stress distribution corresponding to the 
bending of a uniform bar of width a, for which the bending couple is equal 
to 0-5ka?. Hence, the constraint factor 2M/ka? is certainly not greater 
than 1-38 for all notches. This estimate will now be improved by taking 
into account the geometry of the notch and by considering stress and 


velocity fields in some detail. 


(a) Dee p note h with a circular root 


Suppose first that one side of the bar is plane and that on the opposite 
side there is a notch whose root is a circular arc, centre O and radius r 
Fig. 3(a)). To fix ideas it will be assumed that the bar is bent to open the 
notch so that the stress, o,,, across the minimum section is tensile near 
the circular root and compressive on the opposite side. The slip lines in the 
field LN L, defined by the stress-free surface LL, are logarithmic spirals. 
On the surface LL the mean compressive stress p = —k and, according to 
Hencky’s stress relations along the slip lines, p decreases as we travel 
inwards along any slip line. Hence, a,, rises steadily from the value 2k at 
the curved surface. Near the plane side the slip lines are straight and there 
isa uniform compressive stress — 2k parallel to the surface. 

These two slip-line domains may be extended to meet at a neutral point A 
onthe minimum section (Fig. 3 (a@)). N isa point of stress singularity, p being 
discontinuous; o,, changes abruptly from a tension to a compression across 
Y. The position of NV depends on the geometrical parameter r/a and is 
determined by the condition that there is no resultant longitudinal tension 
in the bar. Since V is a neutral point and not a neutral axis, the mean 
compressive stresses on either side are independent. We shall see later, 
however, that if their difference, Ap, is greater than kz the yield condition 
is exceeded in the rigid regions 

It was shown in the preceding section that there cannot be a velocity 
discontinuity along any slip line through a point such as N. Therefore the 
only compatible velocity solution consists of a rotation of the rigid regions 


R. Hill, private communication. 
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about NV. By comparing the slip-line and transformed-velocity fields it js 
easily demonstrated that the rate of plastic work is everywhere positive. 
Another solution satisfying stress and velocity boundary conditions can 
be constructed in the following manner. The two slip-line domains defined 
by the stress-free surfaces are connected by two curved slip lines PQ 











Fic. 3. Slip-line fields for a single notch with a circular root. 


(Fig. 3(b)). The material inside PPQQ is assumed to be rigidly constrained. 
Outside this pivot the rigid ends of the bar rotate by sliding over PQ which, 
therefore, must be a circular arc. The angle which the slip lines PQ and 
PL make with each other must be at least $7 to avoid exceeding the yield 
condition in the angle of the rigid region at P. It cannot, however, be greate! 
than $7 since this would lead to velocity discontinuities along LP and MP 
involving negative plastic work; the situation is similar, though not identi- 
cal, to that in Fig. 2. The same arguments apply at Q@. Hence M PQR must 
be a continuous slip line. It is easily verified that the rate of plastic work 
is then positive throughout the field. 

Let us denote the angle LOM by 2%, the radius of the are PQ by R, and 


its angular span by 2A. The values of , A, and R/a are sufficient to fix the 


shape of the field, and hence the distribution of stress, for a given value ofr 4. 
Following the variation of p along M PQR, we find, by Hencky’s theorem 


that 2 = lta 


This relation, together with two equilibrium equations expressing the con 


dition that the load applied to each half of the bar is a pure couple, enable 


the three unknowns to be calculated. The couple M is most convenient) 
found by taking moments about the centre of curvature of the are PQ. 
It follows from (3) that the couples calculated from both the above slip 
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line fields are upper bounds. The relations between the constraint factor 
and the ratio a/(a+-r) are shown in Fig. 4. The constraint factor obtained 
from the first field rises steadily from the value unity as a/r increases. The 
second field is only valid for A > }n, since for lower values R/a becomes 
negative. It is identical with the first when R/a 0, corresponding to 
ai(a-+-?r) 0-64; then A 


;7, and the jump Ap in mean compressive stress 
at Niskzw. For greater values of a/r the constraint factor obtained from the 
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Fic. 4. Constraint factors for circular notches. 


second field is the smaller, and reaches a maximum value of about 1-261 
corresponding to vanishingly small root radius (r/a—> 0). It follows from 
the converse of (3) that, in this upper range of a/r, the first incomplete 
solution certainly violates the yield condition at some part in its rigid zone 
which coincides with a deforming region of the second slip-line field; 
calculations show that such coincidence occurs only along a segment of the 
are PQ. It isin this range of a/r that Ap > kz, and one would expect over- 
stressing in the rigid radial singularity about NV, since the Hencky stress 
relations prescribe a change in p of only kz along a slip-line are of angular 
span $7 

Slip-line fields may be constructed in a similar manner for a bar sym- 


metrically notched on opposite sides. Fig 5 (a) shows the first type of field 
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and Fig. 5 (b) the second, both of which are symmetrical about the perpen- | be at 
dicular bisector of the minimum section. The corresponding constraint | strain 
factors for different values of a/(a--r) are shown in Fig. 4. The second type relati 
restri 
.@) 
rf. 
ai 
 _— 
| 
a 
a 
(a) (b) 
Fic. 5. Slip-line fields for a pair of notches with circular roots. 
of field again gives the better upper bound when it is valid (i.e. for A > }7 The 
and a/(a+r) > 0-398). The constraint factor has a maximum value of thos 
1-380 when r/a—> 0 and the field consists only of the two circular slip lines N 
PQ; the mode of deformation is the same as that envisaged in order to | 9” 
obtain 1-38 as the upper bound to the constraint factor for all notches. the | 
The 
(b) Deep we dge -shaped notch 
The methods used to construct slip-line fields for notches with circular 
roots are clearly applicable to any shape of notch. In the vicinity of a 
wedge-shaped notch the slip lines defined by the stress-free surface, OL. Whe 
are straight and meet the surface at 45° (Fig. 6). The singularity at 0 couy 
enables the field to be continued round O as far as necessary to form the bett 
region OM P, consisting of straight lines and concentric circular ares. It T 
may now be extended across the minimum section to meet the field defined tiny 
by the opposite stress-free surface at a neutral point V; POP is a right R's 
angle and the region POPN is a region of uniform tension. Alternatively, not 
the domain LM PO may be connected to the opposite domain by a circular slip 
slip line joined smoothly to OP. The angle POP, and hence 2A also, must doe 
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at least $7 to avoid exceeding the yield condition in the rigidly con- 
strained region below POP. It follows, therefore, from the Hencky stress 
relations that for this second field to be valid the angle of the notch, 28, is 


estricted as follows: 


B < 57-3°, single notch; 
ye (7) 
8 73-6°, pair of notches. 








(a) 


Fic. 6. Slip-line fields for a single wedge-shaped notch. 


The fields constructed in this way for a single notch are shown in Fig. 6; 
those for a pair of opposed notches are similar. 

The calculation of constraint factors for different values of B (Fig. 7) is 
similar to that for circular notches, except that the extent of the field at 
the surface is measured by the length s of OM instead of by an angle w. 
rhe yield-point couples obtained from the first type of field are, 


Vl bka*| 1+-(7—28)/(4+-27—28)], single notch; (8) 
MU = tka*{1+-4n—6], pair of notches. 
When A = }7, R/s = 1 and the two types of field give the same yield-point 


) 


ouple; for the lower ranges of 8, defined by (7), the second type gives the 
better upper bound. 
The constraint factor increases with decreasing 8, and this increase con- 


tinues with the second type of field until a value of 8 is reached for which 


Ris oo. This occurs when tan 2A 2/(2A—1), B = 3-2° for a single 
notch; tanA = 2A, 8 = 30-0° for a pair of notches. Beyond this point the 


slip-line field remains the same (s (0) and an increase in sharpness of notch 


does not increase the constraint factor. Note that the maximum constraint 
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factors (1-261, single notch; 1-38, pair of notches) are the same as thoy: 
obtained for circular notches with vanishingly small root radius. 
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Fic. 7. Constraint factors for wedge-shaped notches. 


(c) Deformation with a deep notch 

The dimensions and constraint factors which have been computed for all 
the above fields are given in Table I. The manner of yielding of a bar with 
a single deep notch can be checked experimentally by examining, after a 
small plastic deformation, the shape of the surface which was originall; 
plane. If the bar yields in the manner envisaged in the first type of field 
(Figs. 3 (a) and 6 (a)) the surface region SS should remain plane with a kink 
at each end marking the beginning of the undeformed regions. If, however 
the deformation is in accordance with the second type of field (Figs. 3() 
and 6 (b)), there should be kinks at S and R separating.segments of surface 
which remain plane except close to R. The angular rotation of the tw: 
regions RS is half that of the rigid ends of the bar. The velocity normal 
to the bar of the region RR is constant aleng its length and greater that 
the normal velocity of RS at R owing to the velocity discontinuity through 
R. Hence, according to whether the bar is bent to open or close the notch, 
either a bulge or hollow should develop at RR, rounding off smoothly into 
the region RS close to R. It is to be expected that in a real material this 
kind of surface deformation will only be observed if the material has a sharp 
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yield-point. The values of SS/a and RR/a for different shaped notches are 
ven in Table I and an account of experimental observations is given in 


the Appendix. 


TABLE | 


Constraint factors and dimensions of slip-line fields for pure bending of 


notched bars 


CIR AR NOTCH | WeEbDGE-SHAPED NOTCH 
Ria |2Mjkat| SSja | RRja| 8B sa | Ria |\2M/ka*| SS/a | RR/a 
1°025 So°2 326 1°079 1°079 
+4 I 54 74°5 311 I'lig I'l19 
N 8 y 7 I 8 1-095 68:8 *298 1°156 1-156 
1-135 | 1°135 63-0 -286 te % I‘1g! I-1gt {1-222 
6°7 »9 I'l 5 "92 7°3 275 1°222 s-3e8 | 4 
e530 ahs 57°3 73 |\-275 L +444 
198 1-168 1°136 -706 49°3 219 “314 1°239 I°I57 "432 
7 5 I°217 r'o8s5 "523 37°3 *147 *353 1*253 1-088 “418 
59°5 359 I*2600 I‘o2I 401 17°32 952 "354 1°2601 1030 "404 
32 Oo *389 1*2606| 1-021 *401 
I*4 I .) S50 *354 1-087 
8 800 354 1°174 
\ 76. +262 
I 75°0 354 Ie | 7°29 
73°¢ "35 1°28 
739 354 L-3545 4 
I°4 204 I 63°6 *234 456 1°345 
$ 347 I*25!I 53°60 “144 “510 I°370 
$59 14 43°60 074 535 1°375 
30°0 o "544 1°380 








d) Shallow notches 

The region of plastic deformation is confined to the neck of the bar only 
if the notch is sufficiently deep; if w is the width of the bar, w/a must be 
greater than some critical value p,, say. A lower limit to p, is given by the 
condition that: (i) A deep notch must extend at least as far as its slip-line 
field at the notch surface , 

Another lower limit can be deduced from an upper bound to the constraint 
factor which takes into account the size of the notch though not its shape. 
The bound is obtained from a deformation consisting of two zones, one of 
uniform tension and the other of uniform compression parallel to the axis 
of the bar, enclosed by two intersecting straight lines inclined at 45° to this 
axis. If A is the total area cut away by notches, and S is the total area 
within the criss-cross lines including A, (2) gives an upper bound of 
28S—A)/a2 to the constraint factor. The smallest value of this is 
w°®—2A4) a2, when S has a minimum value of $w? and the lines cross in the 
centre of the bar. This is valid for any number and shapes of notches lying 
within the deforming zones. It is best suited to shallow notches, and leads 
to the following lower limit to p,: (ii) If (w?—2A)/a? is less than the con- 


straint jac tor fon a deep notch o} similar shape 7 then wia is less than Pe- 
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Calculations of these two lower limits have been made for wedge-shaped 
notches, on the assumption that the slip-line fields given above for deep 
notches of this shape are the true solutions. The results for a single notch 
are shown in Fig. 8. The better value of the bound for values of B < 29 


is that given by (ii), but the maximum value 1-21 is given by (i) when 











LOWER LIMIT TO Pc 





| } | | 


30 60 70 80 90 





O id 20 30 40 2 
_ B 


Fic. 8. Lower bounds to p, for a single wedge-shaped notch. 


> 


8 ~ 57°. The results for a pair of equal and opposite notches are similat 
the condition (ii) being best for 8 < 46° and its maximum value obtained 
from (i) being about 1-30 when 8B ~ 68°. The actual critical depth might 
be expected to decrease steadily as 8 increases, but this is not certain and 


can only be decided by accurate solutions for shallow notches. 


Discussion 

The experiments described in the Appendix tie up well with the theory, 
and it is likely that the slip-line fields presented in this paper are 
the true solutions to the problems. From calculations of the averag 
shear stress over oblique sections of the rigid ‘pivot’ bounded by circular 
ares, which occurs in the second type of field for a pair of sharp wedge- 
shaped notches, it appears that the yield stress is unlikely to be locally 
exceeded. However, the yield-point loads calculated from them must be 
regarded as upper bounds unless statically permissible stress fields can be 
constructed in the rigid regions. 
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The same type of slip-line field may be used to estimate the yield-point 
loads for a variety of combinations of tensile force and bending moment 
acting across the minimum section of the bar. 

Constraint factors for notched bars pulled in tension have been calculated 
by Hill (4) and these are much higher than the estimates obtained in this 
paper for similar bars bent by pure couples. For example, the values ob- 
tained with a pair of parallel-sided notches, each with vanishingly small 
root radius (slit-notch), are 2-571 and 1-380 respectively. High tensile 
stresses, however, may occur in the neck with both systems of loading. 
The maximum tensile stress in the region of deformation in a bar with a 
single slit-notch bent open by pure couples occurs at the root of the notch 
ind its value is 2k x 2-043. An even higher tensile stress may possibly occur 

n the rigidly constrained pivot. For comparison, a bar with a pair of slit- 
notches which is pulled in tension has a maximum tensile stress equal to 


};< 2-571 across its minimum section. 
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APPENDIX 
Comparison with experiment} 
Constraint factors 

Wide notched bars of different material were subjected to four-point loading, so 
at over the middle portion containing the notch the applied load was a pure couple, 
nd the bending tended to close the notch. The three materials used were (a) O.F.H.C, 
pper, cold-rolled 50 per cent. from the annealed state; (b) annealed O.F.H.C. 
pper; (c) annealed 18/8 stainless steel. The bars were notched on one side only, 
e other side being a plane polished surface. The notches were wedge-shaped and 
each material two semi-angles, 8 30° and B 70°, were used. The dimensions 
the bars are given in Table II. As the load was increased, a measure of the 
listortion due to bending was given by the deflexion of a fixed dial gauge, whose 
nt rested against the plane side of the bar opposite the root of the notch. The 
bars were bent far enough to locate the sharp bend in the load-distortion curve. 
yield-point load was taken as the ordinate at the point of intersection of the 
ingents to the straight portions of the load-distortion curve on either side of the 

sharp bend. A total deflexion of about 4° proved sufficient in all cases. 
The yield stress in uniaxial tension, Y, was determined for each material by 
rdinary tensile tests. The maximum shear stress in pure shear at yield, k, was 
taken to be Y/3, in accordance with the Mises yield criterion, which is believed to 


ea good approximation for copper and alloy steels such as stainless steel. 


+ The experiments were carried out by Dr. Hundy of the Metal Working Laboratory, 
B.LS.R.A., Sheffield. 
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The constraint factors derived from these experiments are compared with the 
theoretical values in Table IT. The work-hardened copper, which corresponded most 
closely to the theoretical ideal material, gave constraint factors slightly lower thay 
the theoretical values; this may be because the bars were of finite width and th 


deformation only approximately plane strain. The extrapolation for Y with bot} 
TABLE I] 


Experimental results—Pure bending of bars with a single wedge-shaped notch 
I ) { { } 


Experiments Theor, 
Thickness 
Width at notch | Semi-angle Ratio of Rat 
of bars 1 of notch k Y/<3 | Constraint | constraint | Constraint | constra 
Meta (in (in p° (tons/sq. in. factors factors factors factor 
50%, work-hardened 3"O *293 70 I*I3 1°145 
12°27 1085 1*095 
O.F.H.«( copper "Oo *305 30 I*23 I°255 
Annealed 2°0 155 70 Ils 1148 
3°26 1'07 I 
O.F.H.C., « 1! 2* 155 30 1°23 1*255 
Annealed "0 "299 70 1°26 1-148 
13°3 I'l 1095 
18/8 stainless steel 8) *203 30 1°39 1°255 


the annealed materials was rather uncertain due to the rounding of the tensil 
stress-strain curves. Not too much significance, therefore, should be attached to 
the constraint factors obtained with either of these materials. It is interesting to 
note that the more reliable ratio of the two constraint factors, for the sharp and 
wide notches, agrees closely with theory for each of the three materials. 


(b) Deformation 

Taper sections were made of the plane surface of each bar opposite the notch after 
bending, in order to reveal the surface deformation. This was well defined on the 
work-hardened copper, and the tracings of these taper sections (Fig. 9) show how 
close was the agreement with theory. The surface deformation of the annealed 


materials was similar but more diffuse, as one would expect. 


‘01’ 





1° 
B=10° 
<34” 





Fic. 9. Taper sections, work-hardened copper 
(theoretical figures in brackets). 

Hardness tests on cross-sections of the annealed copper specimens after bending 
revealed that the maximum deformation took place at the roots of the notches, as 
one would expect. The hardness readings were too scattered to indicate clearly the 
whole region of deformation, but there was rough agreement with the theory. 

The deformation inside the bars was best illustrated by etching bent mild steel 
specimens in Wazeau’s reagent to reveal the plastically deformed regions. The 


photographs of the etched cross-sections shown in Fig. 10 clearly demonstrate that 
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e type of deformation was the same as that envisaged in the theory (see Figs. 6 (a) 


. 10. Etched cross-sections of mild steel bars showing regions of plastic 
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SUMMARY 

A circular compressible cylinder is given a finite extension increasing its length 
in the ratio A: 1 and then a twist % per unit length so that plane sections perpendicular 
to the axis remain plane sections. The strain-energy function has a completely general 
form in terms of the three invariants of strain and the solution is developed as a 
power series in ys. The ratio of the radii of the strained and unstrained cylinder and 
the longitudinal forces required to maintain the extension and twist are found ip 
series as far as the terms in ys”, whilst the axial couple necessary to produce the twist 
is given as far as the term in Ys. The special case in which the longitudinal force js 
zero is also considered and the resulting extension ratio A is evaluated in this case 
for a given 


1. Introduction 

A NuMBER of problems of finite deformations of isotropic incompressible 
bodies have been solved by Rivlin (1, 2, 3), Green and Shield (4, 5), and 
Green, Rivlin, and Shield (6), using a completely general form for the strain- 
energy function. When the body is compressible, problems have mostly 
been solved by assuming special forms for the strain-energy function and 
progress has, perhaps, been limited because of uncertainty about the precise 
form to be chosen for the strain-energy. The comparatively simple problem 
of a general homogeneous deformation has been solved by Rivlin (1) for 
compressible bodies using a general strain-energy function. Some progress 
on general lines has also been made by Green, Rivlin, and Shield (6) for 
general compressible materials, although solutions of problems were re- 
stricted to the incompressible case. 

It appears to be of theoretical interest to obtain the solution of other 
problems for compressible bodies, making no assumptions about the strain- 
energy function. Such solutions will have practical applications imme- 
diately the special form of the strain-energy function is chosen. In this 
paper the finite extension and torsion of a circular compressible elastic 
isotropic cylinder is considered, making no special assumptions about the 
strain-energy function. Two cases are discussed: 

(i) The cylinder is given a finite extension along its length, of extension 

ratio A, and a finite twist in which each section perpendicular to the 


[Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 2 (1953)] 
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xis remains a plane section perpendicular to the axis after torsion, 
ind rotates through an angle #z proportional to its distance z from one 
end of the cylinder in the strained state after finite extension. 

ii) The cylinder is given a finite twist about its axis, but the resulting 
extension is determined so that the resultant longitudinal force on 
the ends of the cylinder is zero. 

In both cases we assume the displacement and we find the forces and couples 
which are required to maintain the displacement. Problem (ii) has been 
studied by Kappus (7), Riz (8), and Seth (9) using special stress-strain 
relations, but Seth’s results cannot be deduced from the general theory of 
the present paper since (as has been pointed out by Rivlin) his stress—strain 
relations are not derivable from a strain-energy function which is a function 
of the strain invariants only. The solution is given in the form of a power 
series in /, but only terms up to the order of ? are evaluated here. To this 
order of approximation the series solution could have been derived by using 
Murnaghan’s (10) second-order approximation for the strain-energy func- 
tion.t Further terms in our series solution can be obtained if required by 
using the general approach indicated in this paper. 

\ brief summary of the notation (see 5, 11 for full details) will now be 
given. A curvilinear system of coordinates 6; is defined in the body which 
moves with it as it is deformed. In this the Latin suffixes take the values 
1,2. 3. The covariant and contravariant metric tensors for the coordinate 
system 6; in the unstrained body are denoted by g;;, and g“ respectively, 
and for the coordinate system in the strained body the corresponding metric 
tensors are G.,.and GG. whilst 

9=\9unl, G= |Gyl. (1.1) 
The strain invariants are 
G 
L= 9G, i, = Leg, @", I, . (1.2) 
g 
and the contravariant stress tensor (per unit area of the strained body) 


referred to the coordinates @; in the strained body is expressed in the form 


T g*eOo+ BP G*p, (1.3) 
ew OW OW 
wher ® = 21; . P28, p= 2h; (1.4) 
" Ol > al Fel 
— = 3 
Rik g* I, —grg*G,.. (1.5) 
With this notation the equations of equilibrium become 
aik . 
Tr“; = 9, (1.6) 
Note add proof. The writers have now seen a book by F. D. Murnaghan on 
Finite Deformation of an Elastic Solid’ in which he gives a partial solution of case (ii). 


R 


































242 





A. E. GREEN AND E. W. WILKES 


where 7**||,; is the covariant derivative of 7‘* with respect to the strained 
body. The boundary conditions for boundaries at which surface tractiong 
are prescribed are oa 
riky pr. (1.7) 
where n; are the components of the unit vector normal to the bounding 
surface, and P* are the components of applied surface traction. The physical 
components of the stress tensor, along coordinate curves in the strained 
> ano atk Y Yii\s 
body , are T*(Gu G"'):, 


2. General solution 

In the present problem a point in the strained cylinder is specified by 
means of rectangular cartesian coordinates (y,, Y>. ¥3) in which the y,-axis 
is taken along the axis of the cylinder and the origin is at one end. In the 
strained body, cylindrical polar coordinates (7,0,z) are taken to coincide 
with the curvilinear system (6,, 4, #;), and in the strained state the cylinder 
occupies the region 0 < z <1,r <a. It follows that 


Y, = reos8, Yo = rsin8@, Me == 2. (2.1) 


From the assumptions about the deformation, the rectangular coordinates 
(1, ¥g, %3) of a point in the undeformed body are given by 

x, = 71Q cos(O—wz), %, = rQsin(d—wyz), Xe i’ (2.2) 
where the coordinate systems x; and y; coincide, and Q is a function of r 
and 5? to be determined. 


The metric tensors G;,, G*, and q;,, g'* are thus 


Gi, 1 0 O\, Gik 1 0 OV, G = r*. 
0 7 0 0 r-2 0 (2.3 
o> 0 I 0 0 7 
Vix (Q+7rQ,)? 0 0 
0 r2Q? reQry 
0 — rb \-2 1 72Q2y2, 
(2.4) 
git (Q+7rQ,)? 0 0 \, g=A*r’Q?(Q+7Q,) 
0 r 20) 2 Ay? rus 
0 AY ? 


The Christoffel symbols of the second kind for the strained body are all zero 


except that : a 
m,=-—r, IT, 
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r-1, 
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The three strain invariants are 

I, 9G. (Q rQ.) 21Q 2] rap? 1 r2 
975 G1, [((O+7rQ,)?+ Q?+A*+ 1? Q7y? |, | 


i. 
5 a (2.5) 
G 
I, 2Q-2(Q+rQ,)? | 
g 
nd the tensor B“* becomes 
Bu (Q } rQ.) 2(Q 2 i r2 | A? ry") 
R22 (Q rQ.) 2(7 20) 2.) d*ys?) | rr 20) 2 | 
L (2.6) 
B33 — 02(0+-7rQ,)?+A°Q-? | 
Bs b(Y+rQ,)-?, BY Bis 0 | 


The components of stress may now be derived from (1.3), and if these are 
ubstituted in the equations of equilibrium (1.6) a non-linear differential 
quation is obtained for Q. It seems unlikely that a solution of this equation 
in be found in general form, but it is reasonable to assume that Q may be 
expanded as a power series in the parameter %. If this assumption is made 
tis theoretically possible to obtain the terms in this expansion successively, 
using also the boundary conditions on the surface of the cylinder. Attention 
s here confined to the first two terms of the expansion, but, apart from the 
complexity of the algebra, further terms may be obtained if required, with- 
out making any special assumptions about the strain-energy function W. 
Since Y is unaltered if y is changed to —, the expansion of Q is assumed 
0 be ptt (r)+..., 

“ 
where is a constant and f(r) a function of r, both to be determined. 


If this value for Q is substituted in the previous expressions, so that only 


the second powers of ys are retained, the following expansions result: 





I, = X24 2y2+ f2(d2v2— 4y3f—Qy3rf’)+... 
I, = 2A? + t+ [2p ?r?— 43 (A? + pw?) f— 2y3(A2+ *)rf’]+... (3 (2.7) 
I, = N24 2d2p.5?(2f-+-rf’)+... } 
gl = p®?—2y3p?(f+rf’)+... 
92? = p22 (2 Qy8r2f) +... £ (2.8) 
g® = 2, g™* = rf, g* = g* = 0 | 


pul pP(A® + ?) + 2p[ A2r? — 2 (A? + Qu?) f—2Qy(A2 pe )rf’| +... 
BPP Pr 2(A2 + 2) + Pr 2h? A2v2— Qu (A2+ Qu?) f—Qy3rf’]+... 
BS — 2)2yu2— 2)2.3y2(2f+-rf’)- 

B82 — 21,24 22 u3f3( f+-rf’)+..., Be — BY — 0. 


It should be noted that g** is exact and B* is correct to the order wb. 
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Now from (1.4), since W is a function of J,, J, J, it must follow that 
®, Y’, p, are functions of J,, I,, J, which can be assumed to be expanded as 
a series in ?, thus 


® = H+yP{A (Ar? —4y3f—2y rf’) + 
+ FLA p27? — 43(A2+ py?) f—2u3(A2+ p2)rf’] | 
2A w?(H— $15 *H)(2f+97f")} 4 
Y = K+yp*{F (Ar? — 4y8f—2ySrf’)+ | 
+ BLA p?r? —45(2+ pu) f—2p3 (0+ 2 )rf'] a 
- 20*u5(D—4I3 1K)(2f+rf’)}4 
p = w+yr{l, E(*r?—4y3f—2u3rf’)+- 
+ I; Dia? u?r?— 43 (A? +?) f— 23 (02+ pw? rf] 
— 25 (I, C+-415 'w)(2f+7f’)}+... 


where 
H Ee K 2g w ae | 
1 = 2]; oF B = 2]; rae C = 21; Oe t (2.11) 
| 
D = 21; Pa B= 2; ar 5? PF = 21; FFE | 


all being evaluated at 4% = 0. Since, from (2.7), when % = 0, J, Jy, J, are 
all constants depending only on A and yp, it follows that A, B, C, D, E, F,H, 
K, m are all constants. Finally the above expressions are substituted in 
(1.3) to obtain the expansions for 7“ in terms of nine constants arising from 
the strain-energy function W. 
We find 
gil pw? +-?(A? ; p*)K bay ta? (A A272 4 B'f+ C'rf’) a 
722 r 2 2H + 2(A2 +p?) K 4 w | } | 
| 
»—2,),2 mT 2y2_1 rf i. DF MN \r fk’ 11 | 
+7 ub [(A H)A ? B'f | (B ( \rf vee | 2 12) 
733 — H+ 22? K + w+ y?(F'A*r?+ G'f+4G4'rf’)+... 
732 — \*yb(H+-p?2K) +A D’A2r? + E'f+ (4.4 —p3K)rf"|+... | 


7i2 7is 0 
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where 

A’ = p2A+p4(?+p2)B+A2u8 D+ Aut E+ p2(A2+ 22) F 4+ 2K } 
—B’ = 4u5A+ 4y5(A2+ p?)? B+ 4A4u9C + 8A2n7(A2 + ?) D+ 

L 8\2u7 B+ 8u5(A2+ 2) F+ 2u5K 4+ Que 
0! = 4B’—p3H—)y3K 
D' = A+p4B+2y?F +, (2.13) 
BE’ = 2uH+4y3K—4yA’ 
F’ = A+ 2024 B+ Aue D+ Aut 322k 
G’ = 402u3.A + 8)2u5(A2+ 2) B+ 4A4u9C + 402 7(3A2+ 2) D+ 
+ 402.5(A2+ pw?) B+ 402 u3(A2+ 32) F— 22H + 2a | 





The boundary conditions (1.7) are now applied at the curved surface 
r=a. These give 74 = 0, since n; = (1,0,0), and 7!* = 713 = 0. Hence, 
on equating to zero the powers of y*, the following relations are obtained: 


w'K+p?(H+A°K)+ a7 = 0, (2.14) 
A'\a?+- B'f (a)+ C'af'(a) = 0. (2.15) 


Equation (2.14) determines » when A is given, and in virtue of this relation 
71 and 722 take the values 

rll — Y2(A'N224 B+ O'rf’) 4+... \ 

722 — 24° (A’4+ H)d2r2+ Bf + (B’—O'rf']+... 0 

The equations of equilibrium (1.6), 7||; = 0, are already satisfied for 

k = 2, 3, since 7!* = 78 — 0, and the remaining equation for k = 1 becomes 


U4 [2,4 Ph, 722 = 0, 
, dr rit 9s 
that is + —__.__ rr*2 — 0 
dr r 

or C'rf" +3C'f'+A*r(2A’—H) = 0. (2.17) 
The solution of this equation is 

» A272 ’ , ° 

f(r) apy (24 —H)+ar-?+8, (2.18) 


where « and f are arbitrary constants. In this problem the constant « must 
be zero to give a finite solution on the axis, and the constant f is determined 
from equation (2.15), giving 
\?a*{ 2.4’ B’—4A'C’—H(B’'+2C’ 
~— SES eee (2.19) 
8B’C 
\a*| (H +-A?K)(2A'+H)p3— H B’| 


4B'C" 





R 
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The solution (2.18) also solves the problem of the extension and torsion of 
circular cylindrical tube deformed in the same manner as the cylinde; 
described in the present problem. In this case a and f are determined by 
two boundary conditions at inner and outer curved surfaces corresponding 
to (2.15). 


3. Given finite extension and twist 


The boundary conditions (1.7) give, for the surface traction components 
at the ends of the cylinder, the values 


Pi — Q. p2 — 722. Ps — 733. 


since the unit vector in the direction of the normal is n, (0,0,1). The 
corresponding physical surface traction components are (0, P?\ Gy, P?VG,.), 
that is (0, r73?, 73), 

The resultant couple about the axis at the ends of the cylinder necessary 
to maintain equilibrium is given by 


a 


M = 2r | r?(r732) dr. 
0 
When the value of f(r) is substituted from (2.18) in the expression for 7 
from (2.12) and the integration performed the resultant couple becomes 
mA2atgs 
5 | 


MN = H+-p?K)+ 
aAta%y3 


3g p77 EB 24" —H)?—3C' E'(2A'+-H)+8B'C'D'). (3.1) 


The resultant force at the ends of the cylinder necessary to maintain 
equilibrium is obtained from 


a 


N = 2m | rr®™ dr. 
0 
This integral is evaluated, as in the case of the couple, on substitution for 
f(r) and 7*8 and gives 
. ae — TAA ay? cal P 
N = 7a?(?H + 222K + a)-+4 — [2B'F'—G'(2A'+H)]. (3.2) 
‘Zn * 

The couple M in (3.1) is correct to the order ¥®, and the resultant force 
in (3.2) is correct to the order ?, and these expressions are completely 
general for an isotropic compressible body. Further terms in these expan: 
sions in powers of ys may be obtained by a similar process. 
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4, Given twist and zero longitudinal resultant force 
When the cylinder is twisted by a couple alone and is allowed to extend 
so that the resultant force along its axis is zero, then the extension ratios 
\and » can be determined in terms of the twist y%. In the undeformed state 
of the cylinder A = p |, and so it is assumed that 
A 1+A’'s?-+-..., pe 1+ pf? +..., 
where A’, «’ are constants to be determined. If these expressions are now 
substituted in equations (2.14) and (3.2) taking V 0, and account is taken 
of the fact that all the quantities H, K, w, A, B, etc., contain A and yp, it is 
possible to determine A’ and yw’ and hence obtain M from (3.1). It is more 
expedient, however, to revert to the beginning and follow through the work 


by assuming 


A = 14A'pP+..., Q = 14+7f(r)+.... (4.1) 

The quantities H, K, w, A, etc., will now be required when A = 1, p = 1 

in formulae (2.10) and (2.11), that is for 4 = J, = 3, J, = 1, and it will be 

understood henceforward that when the above letters H, K, ete. are used 

they carry this implication. 

It now follows from (2.7) that 

I, 3+?(r?—4f—2rf’+2X')+... 

I, = 3+2(r?—8f—4rf’+4)')+... }, (4.2) 
I, = 1—2p?(2f+rf’—A’)+... | 


whilst from (2.8) and (2.9) we obtain 


gi 1 —2?( f+rf’)+. | 
g? = 24222 2f) | ; 
g°3 = 142)’ — 
ie b+ 2A" q* g™* 0 
Bu — 24 .42(r?—6f—4rf’+22’)-4 
B= 2r-*-+- 1"? (7? — 6f — 2rf’+-2X')+... 

: (4.4) 
B »— 2y?2(2f+-rf’ —2N')-+... 
B? = f—2F( f+rf'—d’)+..., BY — Be=0 


Further, from (2.10) it follows that 
® = O'+A'f?(2A+4F+42E—H)-4... 
Y= ¥'+1'?(2F +4B+2D—K)-+... }, (4.5) 
p = p’+N2(2E+4D+2C+a)+... | 


where @’. Y’ 


p’ are the values of ®, VY’, p given in (2.10) after A = p ] 
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have been substituted. From the above, new expressions for 7 are obtained 
on using equation (1.3). Thus 


gil — 7/N_4@’)'f?+..., (4.6) 
722 — 7'22_1G')'r-2y24 ..,, (4.7) 
733 — 733 _ CN if? ..., (4.8) 
732 — 7/824 (24’1+ H—K)X'f-+..., (4.9) 


where 7’1!, etc., are values taken from (2.12) with A and p set equal to unity, 
The values of A’, B’, etc., are those given in (2.13) after substitution of 
A = p = I, and are given below: : 
A' = A+2B+D+£E+3F4+K 
— B’ = 4A4+16B+4C0+16D+8E+16F+2K+20 
C’ = 4B’—-H-K 
D' = A+B+2F L (4,10) 
E’ = 2H+4K—4A' | 
F’ = A+2B+D+£E+3F | 
—G’ = 44+16B+40+16D+8E£+16F—2H+ 20 | 
From these equations it follows that 
B’ = G’—2H—2K, F’ = A’—K. 
The boundary condition at the surface r = a is 71! = 0. Since this must 
be true for all values of % the following two equations result: 


H+2K+ao = 0, (4.11) 
A'a*+ B'f (a)+C’af’(a)—43G'N"' = 0. (4.12) 
The equation of equilibrium gives 
C'rf"+3C'f'+(2A'—H)r = 0, (4.13) 
with the solution f(r) = - <7; (24'—H) +B. (4.14) 


The constant f is determined from the equation (4.12) and yields 
a*|(H-+K)(2A’+H)—HB’] , @’N’ 
Be aB SOB 
The additional information required in this section for the evaluation 
of A’ is obtained by putting the resultant longitudinal force N equal to zero. 
Thus 


(4.15) 


a 
N = 22 | 73 dr = 2m | r(r'3—C'D'f?) dr 
0 0 
rays? 


y/9 


= [2B F’—@'2'+H)]—natnye( 0" — Fe) — 0, (4.16) 


2B’ 
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in which equations (4.8) and (4.11) have been used. This equation leads to 
the value of A’, namely 
a*{2B' F’—G@'(2A’+H)] 


X’ ER aD cits Be 


2(G’2— 2 B’C’) ei 


The axial couple necessary to hold the twist is then 


M 2a [ 73732 dr 


Tr 4); i 7r¢ 6,3 ’ , ‘ ae Al , ' mune 
weet) 4 eB (2A'—H)?—30' E'(2A' +H) +8B'C'D'}— 











7a5y3| G’ BE’ +-2B'(2A'+H—K)|[2B’ F’—G'(2A’+H)] (4.18) 
8.B'(G’2—2B’C’) : ; 
In addition, the radius a, of the unstrained cylinder is given by 
A all - us*f (a) i —_ 
a**{(2A’+H)C’—F'G’|, — 
a 7 oa 
2(G’2—2B’C’) 
ays? . B'H 
a\1l- F 2K —H—4A'+ B’ —-———_} +... |. 
‘| 4(3.B'+2H-4 5K)| ETE) | 
(4.19) 


It is interesting to note here that this series in general involves y? and 
higher powers of %?. This is contrary to the result obtained by Seth (9) 
whose result shows a change in the radius of the strained cylinder dependent 


only on #* and higher powers. 
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ON DYNAMICAL PROBLEMS OF PLANE STRESS AND 
PLANE STRAIN 


By R. E. D. BISHOP (Ministry of Supply)t 
[Received 8 May 1952] 
SUMMARY 
The theory of potential functions in dynamical elasticity is developed. The 
special cases of plane-stress and plane-strain motion are examined, and it is shown 
that the various well-known plane-strain solutions have counterparts in plane stress, 
These correspond, effectively, to a change in one elastic constant (the medium being 
assumed to be homogeneous and isotropic). The representative problem of sym- 
metrical waves in an elastic plate is treated. 
1. Introduction 
THE equations of motion of an elastic solid whose Lamé constants are \ and 
u may, in the absence of body forces, be written in the form 
a9 
ae ed 
(A+p)V(V.d)+yV"*d = p—, (1) 
ot? 
where d = (u,v, w) is the displacement vector and p is the density, or by 
the identity V(V.d) = VxVxd-+V2d, (2) 
‘ cd = pf 
(A+ 2n)V(V.d)—pV x V xd = p- 


et? © 


Integrating with respect to ¢ and rejecting the arbitrary functions of 


integration, we have d = Wf+0U~xA, (4) 
.s s¢ 

where ¢ = (A+2y)p-' | | (V.d)dtdt, A — (p/p) | | (V xd) dtdt. (5) 
0 0 0 0 

Now V.d = V%Z, 


so that V°¢ represents the dilatation. Again, 
Vxd = VxVxA= V(V.A)—V2A = —V?A, 


so the rotation at any point is given by —V?A/2. It follows from (5) that 
¢ and A must satisfy the wave equations 


ran) i a oo 
et as eA 
— c VA = 


et? et? 


c2 VC , (6) 
where c? = (A+2u)/p, c3 = p/p. It is usually simpler to tackle problems 
of dynamical elasticity by the use of (4) and (6) rather than from (1). 
The present note is concerned with two methods of simplifying the 
analysis. These are effected by postulating the conditions of either plane 
+ Now at the Engineering Laboratory, Cambridge. 


[Quart. Journ. Mech. and Applied Math., Vol. VI, Pt. 2 (1953)] 
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stress or plane strain. Several solutions involving the latter are well known; 
but little attention appears to have been paid to dynamical plane stress, 
apart from a brief treatment by Timoshenko (1). 


2. Plane strain 
In a problem involving rectangular coordinates, if we consider an elastic 

solid whose z-dimension is very large so that w and v are independent of z 

and the displacement w is everywhere zero, the conditions of plane strain 

are met. The vector A may be taken to be (0, 0, %), where # is independent 

of z. If, in addition, f = d(x, y), (4) yields 

Ch . Ob - Ch = bs 


Cx cy cy Cx 


(7) 


and (6) reduce to 
9 > o*dh ‘ ‘ CAs 
c2Vid os co Vin a (8) 
at“ ot 
in which Vj denotes é?/éx*-+-6*/éy?. 
In the solution of a boundary-value problem use is made of the stress- 


strain relations. In this kind of motion these are 


- cu —_~ : ov ~ 
v2 AA }+-2u YY AA-4 2p - : Ze AA 
C2 cy 
; (9) 
— cu Cv — — 
xy n| ; F rz = yz = UV 
cy C2 


where A CUu/CxX—-—cCv cy. 


3. Plane stress 

We now consider the motions of a very thin plate occupying the zy-plane 
under the usual conditions of plane stress. Using rectangular coordinates, 
the stresses 22, 22, and 72 are assumed to vanish not only at the surface but 
throughout the thickness, while the three remaining stresses are supposed 
to be independent of z. From Hooke’s law it is found that 


= war: 2 ~ oa Ov in Ou. ov 
ra N'A‘ Qu Ml yy NA‘ Qu R ry = uf: +3). (10) 
Cx ssi cy CY Cx 
where rN’ = 2pA/(A+ 2p), A’ = du/éx+ dv/oy. 


When these values are substituted into the equations of motion for plane 
stress (written in terms of the non-evanescent stresses), it is found that 


(\’ +m), )a’ uV2(u, v) p 3 (ts, ). (11) 


coy col" 


Mathematically, the equations of motion and stress (excluding that for 














bo 
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2) for plane stress are identical with those for plane strain, the only differ. 
ence being that A’ appears in the former and A in the latter. In effect, then, 
the only ditference is one of elastic constants. Therefore, any dynamical 
plane-strain solution will have a counterpart in plane stress; e.g. ‘Rayleigh’ 
surface waves appear in both theories and are qualitatively the same, 
although their speeds differ. 

It remains to be demonstrated that potential functions, comparable 
with ¢ and #, can be found for plane-stress problems. By analogy with 
equations (7) we may write 

u = mee 9 = =, (12) 
ou Oy oy Ox 
Substituting these expressions into the pair of equations (11), we find that 
they are satisfied provided that « and f satisfy the wave equations 


O20 eB 
efVia = —. c2 V268 — — 13 
oe we 2ViB = oe ~ 


where c7, = (A’+-2y)/p. 


4. Symmetrical waves in a plate 

It is of interest to discuss an illustrative problem in both plane strain and 
plane stress. The one selected is that of symmetrical waves in an infinite 
plate. 

1. Plane strain 

Consider a flat plate of thickness 2b which is infinitely extended in the 
x and z directions and whose mid-plane is associated with the surface y = 0. 
If solutions of the form 


wu = (even function of y).exp[ié(a—ct)] ) 


v = (odd function of y).exp[7é(~—ct) | J — 
are sought, the following functions are found: 
¢ = Acosh(émy)exp|ié(x—ct) | ) (18) 
ys = Bsinh(ény)exp|ié(x—ct) | .? 
where A and B are constants and 
m? = 1—c?/c? | (16) 


9; 


n? = 1—c?/c2 J 
. *“,? _ . 
Applying to ¢ and & the conditions that yy and xy must vanish when 
y = +6, a dispersion equation is found relating the wave number é and 
the phase speed c in terms of A, pw, p, and 6. It is 


tanhgénb —— 4mn (17) 
tanhémb- (1-+-n?)?’ 
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PLANE 


anda dispersion curve for the lowest mode may be drawn using the following 


calculated points (for which Poisson’s ratio is vy = 0-3). 


0°543 11 I*304 1°527 1°Q75 2°197 2°645 
I 1°527 1°434 1°247 I°177 1-080 
3°53 4 5 5°450 x 

0985 62 0°935 0°9274 


» 


The wave-length A is given by 27/€ and when this is large compared with 
he dimension 6, €b tends to zero, so that the hyperbolic terms may be 


replaced by their arguments in (17). The limiting value of the phase velocity 
is then found to be c, 

On the other hand, when £0 is infinitely large, the hyperbolic terms may 
be taken as unity. This leads to the cubic 


c& oct 4/24 16 c 
a c a | 
c§ ch c2 ci c 


9 9 


0, - (18) 


ito note 


which defines the speed of Rayleigh surface waves (= 0-9274c, for v = 0-3). 
Between these two limiting values of c/c,, the dispersion curve suffers a 
fairly abrupt drop. 

These results are all well known, and for further details of the motion 


reference should be made to Lamb’s original solution of the problem (2). 


2. Plane stress 
If. instead of being very large in the z-direction, the plate is now supposed 
to be very small in this direction, it may be seen at once that the theory 
isexactly the same as for plane strain, except that ’ replaces A, c%, is written 
instead of cj, and /? replaces m* where 
12 — 1—c2/c2. (19) 
The dispersion equation is now 
tanh nb 4in (20) 
tanh élb (1+-n?)?’ 
and it may be discussed just as before. As might be expected, the general 
shape of the dispersion curve is similar to that of (17), some calculated 


points (for 1 0-3) being given in the following table. 





0°337 5 o’g16 I°150 I°331 I*505 1°624 1°325 
O 1°583 I°552 1°505 1°443 I°391 I*300 
I 2-615 2-918 3°186 37452 4°093 4°607 x 

‘rss 7 1°025 1°000 0"950 o°950 0°937 o°9162 


The phase speed of very long waves is now found to be the ‘bar velocity’ 
Cy (E p)}*, 


where £ is Young’s modulus; if v = 0-3, cy = 1-6125c,. 
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At the other extreme, very short waves are propagated with a speed that 


is found by replacing the hyperbolic functions of (20) by unity. In this 


way it is found that (18) obtains with c? instead of c?. The relevant root, 


p 
for v 0:3, isc 0-9162c.,. 
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NOTE ON THE DERIVATION BY SIMPLE ALGEBRA 
OF ROUTH’S STABILITY CRITERION FOR THE 
BIOLADRATIC CHARACTERISTIC EQUATION 


By J. MORRIS (Guildford) 


Received 8 May 1952] 


Let the characteristic biquadratic be 
F(m) Am'*+ Bm3+ Cm?+Dm-+ E, (1) 
which, without loss of generality, may be written 
f(m) m'+-bm3+cm?+dm--e: (2) 
and let 6. c, d, e all be real and other than zero. 


In these circumstances f(m) can be factorized into two quadratic factors 


of the form - e 
f(m) (m=+ p,m-+-q,)(m-+ pom-+ qo), (5) 


in which the p’s and q’s are real. Equating coefficients of like powers of 


m in (2) and (3), we have 


Pit Pe2 b, (4) 
Py Pot htd = ¢, (5) 
Py dtPoh = 4, (6) 
7193 >> ©- (7) 
From (4) and (6) we find that 
Pi(q2—th) = —bq, +4. (8) 
Pold2—h) = bq2,—d, (9) 
and thus P; Po(Go—)* bd(q, +42) —@—b"q, 2. (10) 
or, making use of (5) and (7), 
D1 Pol (Go—q,)? + bd} (bed —d* — eb?). (11) 


Hence, referring to (2), if the coefficients b, c, d, e and the expression 
(bed —d?—eb*) are all positive, it follows that 

(i) (p,+p,) and p, p, must each be positive; and hence that p, and p, 
must each be positive. 

(ii) (p,¢2+poq,) and q, g. must each be positive; in which case, since the 
p's are positive, the g’s must likewise be positive. 


[Quart. Journ. Mech. and Applied Math., Vol. VI, Pt.2 (1953)] 
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With these conditions fulfilled real roots of f(m) = 0 must necessarily 


be negative, and conjugate complex roots must necessarily have their real 
parts negative. If, however, (bed —d*—eb?) = 0, then either p, or p, must 
be zero, which presages the presence of a pair of roots which are pure 
imaginaries. 

We have thus derived by comparatively simple algebra the necessary 
and sufficient conditions for the ‘stability’ of the ‘transient motion’ repre- 
sented by a biquadratic ‘characteristic’ equation, which conditions are in 
accord with the classic treatment given by Routh in his Advanced Rigid 
Dynamics (section 287, 5th edition). 
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